Skip to content

Commit

Permalink
update vivax clinical model (#400)
Browse files Browse the repository at this point in the history 
* update vivax clinical model

* update Vivax test

---------

Co-authored-by: Amanda Ross <[email protected]>
  • Loading branch information
@rossaman4 @rossaman
rossaman4 and rossaman authored Jan 9, 2025
1 parent 2847706 commit 4c2308f
Show file tree
Hide file tree
Showing 6 changed files with 17,499 additions and 10,460 deletions.
119 changes: 79 additions & 40 deletions model/Host/WithinHost/WHVivax.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
* Copyright (C) 2005-2021 Swiss Tropical and Public Health Institute
* Copyright (C) 2005-2015 Liverpool School Of Tropical Medicine
* Copyright (C) 2020-2022 University of Basel
*
*
* OpenMalaria is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
Expand Down Expand Up @@ -37,6 +37,8 @@
#include <limits>
#include <cmath>



namespace OM {
namespace WithinHost {

Expand All @@ -55,7 +57,7 @@ struct HypnozoiteReleaseDistribution: private LognormalSampler {

/// Sample the time until next release
SimTime sampleReleaseDelay(LocalRng& rng) const {
double liverStageMaximumDays = 16.0*30.0; // maximum of about 16 months in liver stage
double liverStageMaximumDays = 16.0*30.0; // maximum of about 16 months in liver stage
double delay = numeric_limits<double>::quiet_NaN(); // in days
int count = 0;
int maxcount = 1e3;
Expand Down Expand Up @@ -98,6 +100,7 @@ double pRelapseOneA = numeric_limits<double>::signaling_NaN(),
double pRelapseTwoA = numeric_limits<double>::signaling_NaN(),
pRelapseTwoB = numeric_limits<double>::signaling_NaN();
double pEventIsSevere = numeric_limits<double>::signaling_NaN();
std::string vivaxClinOption = " ";

// Set from healthSystem element:
bool ignoreNoPQ = false;
Expand All @@ -113,6 +116,7 @@ VivaxBrood *sampleBrood = 0;

// ——— individual models ———


// number of hypnozoites per brood:
map<double,int> nHypnozoitesProbMap;
void initNHypnozoites(){
Expand All @@ -128,6 +132,8 @@ void initNHypnozoites(){
nHypnozoitesProbMap[cumP] = n;
}
}


int sampleNHypnozoites(LocalRng& rng){
double x = rng.uniform_01();
// upper_bound finds the first key (cumulative probability) greater than x:
Expand All @@ -136,6 +142,7 @@ int sampleNHypnozoites(LocalRng& rng){
return it->second; // corresponding n
}


// time to hypnozoite release after initial release:
SimTime sampleReleaseDelay(LocalRng& rng){
bool isFirstRelease = true;
Expand All @@ -157,24 +164,25 @@ SimTime sampleReleaseDelay(LocalRng& rng){
VivaxBrood::VivaxBrood( LocalRng& rng, int origin, WHVivax *host ) :
primaryHasStarted( false ),
relapseHasStarted( false ),
relapsebHasStarted( false),
hadEvent( false ),
hadRelapse( false ),
origin(origin)
{
origin(origin)
{
set<SimTime> releases; // used to initialise releaseDates; a set is better to use now but a vector later

// primary blood stage plus hypnozoites (relapses)
releases.insert( sim::nowOrTs0() + latentP );
int numberHypnozoites = sampleNHypnozoites(rng);

for( int i = 0; i < numberHypnozoites; ){
//TODO: why do we have two latent periods (latentP + latentReleaseDays added in sampleReleaseDelay())?
SimTime randomReleaseDelay = sampleReleaseDelay(rng);
SimTime timeToRelease = sim::nowOrTs0() + latentP + randomReleaseDelay;
bool inserted = releases.insert( timeToRelease ).second;
if( inserted ) ++i; // successful
// else: sample clash with an existing release date, so resample
}

// Copy times to the vector, backwards (smallest last):
releaseDates.insert( releaseDates.end(), releases.rbegin(), releases.rend() );

Expand Down Expand Up @@ -202,18 +210,20 @@ void VivaxBrood::checkpoint( ostream& stream ){
bloodStageClearDate & stream;
primaryHasStarted & stream;
relapseHasStarted & stream;
relapsebHasStarted & stream;
hadEvent & stream;
hadRelapse & stream;
origin & stream;
origin & stream;
}
VivaxBrood::VivaxBrood( istream& stream ){
releaseDates & stream;
bloodStageClearDate & stream;
primaryHasStarted & stream;
relapseHasStarted & stream;
relapsebHasStarted & stream;
hadEvent & stream;
hadRelapse & stream;
origin & stream;
origin & stream;
}


Expand All @@ -238,10 +248,13 @@ VivaxBrood::UpdResult VivaxBrood::update(LocalRng& rng){
#endif

// an existing or recently terminated blood stage from the same brood
// protects against a newly released Hypnozoite
//NOTE: this is an immunity effect: should there be no immunity when a blood stage first emerges?
if( bloodStageClearDate + bloodStageProtectionLatency >= sim::ts0() ) continue;

// protects against a newly released hypnozoite for relapse classification B
if( (vivaxClinOption=="B1j" || vivaxClinOption=="B2j") && (bloodStageClearDate + bloodStageProtectionLatency >= sim::ts0()) ) continue;

if( !relapsebHasStarted && relapseHasStarted ){
relapsebHasStarted = true;
result.newRelapsebBS = true;
}
if( !relapseHasStarted && primaryHasStarted ){
relapseHasStarted = true;
result.newRelapseBS = true;
Expand Down Expand Up @@ -289,7 +302,7 @@ void VivaxBrood::treatmentLS(){

WHVivax::WHVivax( LocalRng& rng, double comorbidityFactor ) :
cumPrimInf(0),
treatExpiryLiver(0), treatExpiryBlood(0),
treatExpiryLiver(0), treatExpiryBlood(0),
pEvent( numeric_limits<double>::quiet_NaN() ),
pFirstRelapseEvent( numeric_limits<double>::quiet_NaN() ),
pSevere( 0.0 )
Expand All @@ -315,9 +328,12 @@ WHVivax::~WHVivax(){
#endif
}



double WHVivax::probTransmissionToMosquito(vector<double> &probTransGenotype_i, vector<double> &probTransGenotype_l)const{
assert( WithinHost::Genotypes::N() == 1 );
for(auto inf = infections.begin(); inf != infections.end(); ++inf)
for(auto inf = infections.begin();
inf != infections.end(); ++inf)
{
if( inf->isPatent() ){
// we have gametocytes from at least one brood
Expand All @@ -327,6 +343,7 @@ double WHVivax::probTransmissionToMosquito(vector<double> &probTransGenotype_i,
return 0; // no gametocytes
}


bool WHVivax::summarize(Host::Human& human) const{
if( infections.size() == 0 ) return false; // no infections: not patent, nothing to report
mon::reportStatMHI( mon::MHR_INFECTED_HOSTS, human, 1 );
Expand Down Expand Up @@ -360,51 +377,63 @@ void WHVivax::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i, int &nN

for( int i = 0; i < nNewInfs_l; ++i )
infections.push_back( VivaxBrood( rng, InfectionOrigin::Indigenous, this ) );



// update infections
// NOTE: currently no BSV model
morbidity = Pathogenesis::NONE;
uint32_t oldCumInf = cumPrimInf;
bool treatmentLiver = treatExpiryLiver > sim::ts0();
bool treatmentBlood = treatExpiryBlood > sim::ts0();
double oldpEvent = ( std::isnan(pEvent))? 1.0 : pEvent;
// always use the first relapse probability for following relapses as a factor
double oldpRelapseEvent = ( std::isnan(pFirstRelapseEvent))? 1.0 : pFirstRelapseEvent;
double matImmClin = 1 - (0.90 * exp(-2.53*ageInYears));
//double matImmClin = 1 - (1 * exp(-(0.639*8*ageInYears)/2.53));
auto inf = infections.begin();
while( inf != infections.end() ){
if( treatmentLiver ) inf->treatmentLS();
if( treatmentBlood ) inf->treatmentBS(); // clearnace due to treatment; no protection against reemergence
if( treatmentBlood ) inf->treatmentBS(); // clearance due to treatment; no protection against reemergence
VivaxBrood::UpdResult result = inf->update(rng);
if( result.newPrimaryBS ) cumPrimInf += 1;

if( result.newBS ){
// Sample for each new blood stage infection: the chance of some
// clinical event.
// Sample for each new blood stage infection: the chance of some clinical event.
// model variant: no illness from relapses possible unless there was illness from the primary infection

bool clinicalEvent = false;
if( result.newPrimaryBS ){
// Blood stage is primary. oldCumInf wasn't updated yet.
double pPrimaryInfEvent = pPrimaryA * pPrimaryB / (pPrimaryB+oldCumInf);
//double pPrimaryInfEvent = matImmClin * pPrimaryA * pPrimaryB / (pPrimaryB+oldCumInf);
double pPrimaryInfEvent = matImmClin * pPrimaryA * exp(-pPrimaryB * oldCumInf);
clinicalEvent = rng.bernoulli( pPrimaryInfEvent );
inf->setHadEvent( clinicalEvent );
} else if ( result.newRelapseBS ){
// Blood stage is a relapse. oldCumInf wasn't updated yet.
double pFirstRelapseEvent = oldpEvent * (pRelapseOneA * pRelapseOneB / (pRelapseOneB + (oldCumInf-1)));
clinicalEvent = rng.bernoulli( pFirstRelapseEvent );
inf->setHadRelapse( clinicalEvent );
}else{
// Subtract 1 from oldCumInf to not count the current brood in
// the number of cumulative primary infections.
if( inf->hasHadRelapse() ){
double pSecondRelapseEvent = oldpRelapseEvent * (pRelapseTwoA * pRelapseTwoB / (pRelapseTwoB + (oldCumInf-1)));
clinicalEvent = rng.bernoulli( pSecondRelapseEvent );
}else{
// If the primary infection did not cause an event, there
// is 0 chance of a secondary causing an event in our model.
clinicalEvent = false;
}
}


} else if ( result.newRelapseBS ){
// Blood stage is a relapse. oldCumInf wasn't updated yet. Subtract 1 from oldCumInf not
// to count the current brood in the number of cumulative primary infections

double pFirstRelapseEvent = matImmClin * pRelapseOneA * exp(-pRelapseOneB * (oldCumInf-1));
clinicalEvent = rng.bernoulli( pFirstRelapseEvent );
inf->setHadRelapse( clinicalEvent );

} else if ( result.newRelapsebBS ){

if (vivaxClinOption=="A1j" || vivaxClinOption=="A2j"){
double pFirstRelapseEvent = matImmClin * pRelapseOneA * exp(-pRelapseOneB * (oldCumInf-1));
clinicalEvent = rng.bernoulli( pFirstRelapseEvent );
inf->setHadRelapse( clinicalEvent );
}
if (vivaxClinOption=="B1j" || vivaxClinOption=="B2j"){
double pSecondRelapseEvent = matImmClin * pRelapseTwoA * exp(-pRelapseTwoB * (oldCumInf-1));
clinicalEvent = rng.bernoulli( pSecondRelapseEvent );
inf->setHadRelapse( clinicalEvent );
}

} else {
double pSecondRelapseEvent = matImmClin * pRelapseTwoA * exp(-pRelapseTwoB * (oldCumInf-1));
clinicalEvent = rng.bernoulli( pSecondRelapseEvent );
}


if( clinicalEvent ){
pSevere = pSevere + (1.0 - pSevere) * pEventIsSevere;
if( rng.bernoulli( pEventIsSevere ) )
Expand All @@ -418,7 +447,7 @@ void WHVivax::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i, int &nN
else ++inf;
}

//TODO were pEvent and pFirstRelapseEvent meant to get updated?


//NOTE: currently we don't model co-infection or indirect deaths
if( morbidity == Pathogenesis::NONE ){
Expand Down Expand Up @@ -581,7 +610,17 @@ void WHVivax::init( const OM::Parameters& parameters, const scnXml::Model& model
pRelapseTwoA = ce.getPRelapseTwoPlus().getA();
pRelapseTwoB = ce.getPRelapseTwoPlus().getB();
pEventIsSevere = ce.getPEventIsSevere().getValue();
vivaxClinOption = ce.getVivaxClinOption();

// Define accepted values for vivaxClinOption
const std::set<std::string> acceptedOptions = {"A1j", "A2j", "B1j", "B2j"};

// Check if vivaxClinOption is valid
if (acceptedOptions.find(vivaxClinOption) == acceptedOptions.end()) {
throw util::xml_scenario_error("Invalid vivaxClinOption: " + vivaxClinOption +
". Accepted values are: A1j, A2j, B1j, B2j.");
}

initNHypnozoites();
Pathogenesis::PathogenesisModel::init( parameters, model.getClinical(), true );
}
Expand Down
38 changes: 22 additions & 16 deletions model/Host/WithinHost/WHVivax.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
* Copyright (C) 2005-2021 Swiss Tropical and Public Health Institute
* Copyright (C) 2005-2015 Liverpool School Of Tropical Medicine
* Copyright (C) 2020-2022 University of Basel
*
*
* OpenMalaria is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or (at
Expand Down Expand Up @@ -66,8 +66,8 @@ class VivaxBrood{
VivaxBrood( istream& stream );

struct UpdResult{
UpdResult() : newPrimaryBS(false), newRelapseBS(false), newBS(false) {}
bool newPrimaryBS, newRelapseBS, newBS, isFinished;
UpdResult() : newPrimaryBS(false), newRelapseBS(false), newRelapsebBS(false), newBS(false) {}
bool newPrimaryBS, newRelapseBS, newRelapsebBS, newBS, isFinished;
};
/**
* Do per time step update: remove finished blood stage infections and act
Expand Down Expand Up @@ -113,22 +113,23 @@ class VivaxBrood{

// Whether the relapse blood stage infection has started
bool relapseHasStarted;

//Whether the later category of relapses has started
bool relapsebHasStarted;

// Whether any clinical event has been triggered by this brood
bool hadEvent;

// Whether a relapse event has been triggered by this brood
bool hadRelapse;

// Infection origin
int origin;
// Infection origin
int origin;
};


/**
* Implementation of a very basic vivax model.
*
* This is for tropical Vivax (low transmission) and where there is little
* immunity.
* Implementation of a basic vivax model
*/
class WHVivax : public WHInterface {
public:
Expand All @@ -155,19 +156,22 @@ class WHVivax : public WHInterface {
virtual bool summarize(Host::Human& human) const;

virtual void importInfection(LocalRng& rng, int origin);

virtual void update(Host::Human &human, LocalRng& rng, int &nNewInfs_i, int &nNewInfs_l,


virtual void update(Host::Human &human, LocalRng& rng, int &nNewInfs_i, int &nNewInfs_l,
vector<double>& genotype_weights_i, vector<double>& genotype_weights_l, double ageInYears);

virtual bool diagnosticResult( LocalRng& rng, const Diagnostic& diagnostic ) const;

virtual Pathogenesis::StatePair determineMorbidity( Host::Human& human, double ageYears, bool );

virtual void clearImmunity();

virtual InfectionOrigin getInfectionType() const {
return InfectionOrigin::Indigenous;
}



protected:
virtual void treatment( Host::Human& human, TreatmentId treatId );
Expand All @@ -190,7 +194,7 @@ class WHVivax : public WHInterface {

list<VivaxBrood> infections;

/* Is flagged as never getting PQ: this is a heteogeneity factor. Example:
/* Is flagged as never getting PQ: this is a heterogeneity factor. Example:
* Set to zero if everyone can get PQ, 0.5 if females can't get PQ and
* males aren't tested (i.e. all can get it) or (1+p)/2 where p is the
* chance a male being tested and found to be G6PD deficient. */
Expand All @@ -205,10 +209,12 @@ class WHVivax : public WHInterface {
// Either never or expiry time of treatment
SimTime treatExpiryLiver, treatExpiryBlood;

// The probability of a clinical event from a first infection
// The probability of a clinical event from a primary infection
double pEvent;
// The probability of a clinical event from a relapse
// The probability of a clinical event from an early relapse
double pFirstRelapseEvent;
// The probability of a clinical event from a later relapse
double pSecondRelapseEvent;

// Used for reporting; updated by update()
double pSevere;
Expand Down
Loading

0 comments on commit 4c2308f

Please sign in to comment.