Add per genotype vaccine efficacy (#339)

* change default branch in build.sh script

* add Human argument to WHInterface update method

* move nNewNumInfections reporting to after WithinHost update

* add new PEV behaviour and model option switch

* avoid rng call when vaccine is not active

* add genotypes to vaccines

* fix number of new infections

* fix number of new infections in DescriptiveWithinHost (breaks tests)

* fix number of infections to be consistent with the old model

* add genotype to bsv vaccine intervention

* rename PEV_GENOTYPE option into VACCINE_GENOTYPE

* move tbvFactor out of whm.probTransmissionToMosquito

* add genotype option to TBV vaccine intervention

* fix number of infections with VACCINE_GENOTYPE option

* remove debug outputs

* add additional checks in vaccine xml initialization

* add vaccine tests and bump version to 44.0

model/Host/Human.cpp

@@ -128,8 +128,9 @@ void Human::update(Transmission::TransmissionModel& transmission) {
     int nNewInfs = infIncidence->numNewInfections( *this, EIR );
 
     // ageYears1 used when medicating drugs (small effect) and in immunity model (which was parameterised for it)
-    withinHostModel->update(m_rng, nNewInfs, EIR_per_genotype, ageYears1,
-            _vaccine.getFactor(interventions::Vaccine::BSV));
+    withinHostModel->update(*this, m_rng, nNewInfs, EIR_per_genotype, ageYears1);
+
+    infIncidence->reportNumNewInfections(*this, nNewInfs);
 
     // ageYears1 used to get case fatality and sequelae probabilities, determine pathogenesis
     clinicalModel->update( *this, ageYears1, age0 == sim::zero() );

model/Host/InfectionIncidenceModel.cpp

@@ -61,7 +61,7 @@ double EstarInv;
 
 // model options:
 bool opt_neg_bin_mass_action = false, opt_lognormal_mass_action = false,
-        opt_no_pre_erythrocytic = false, opt_any_het = false;
+        opt_no_pre_erythrocytic = false, opt_any_het = false, opt_vaccine_genotype = false;
 
 // ———  variables  ———
 int InfectionIncidenceModel::ctsNewInfections = 0;
@@ -93,6 +93,8 @@ void InfectionIncidenceModel::init ( const Parameters& parameters ) {
 
     opt_no_pre_erythrocytic = util::ModelOptions::option (util::NO_PRE_ERYTHROCYTIC);
     opt_neg_bin_mass_action = util::ModelOptions::option (util::NEGATIVE_BINOMIAL_MASS_ACTION);
+    opt_vaccine_genotype = util::ModelOptions::option (util::VACCINE_GENOTYPE);
+
     if (opt_neg_bin_mass_action) {
         inf_rate_shape_param = (baseline_avail_shape_param+1.0) / (r_square_Gamma*baseline_avail_shape_param - 1.0);
         inf_rate_shape_param=std::max(inf_rate_shape_param, 0.0);
@@ -221,8 +223,11 @@ int InfectionIncidenceModel::numNewInfections (Human& human, double effectiveEIR
     throw TRACED_EXCEPTION (out.str(), util::Error::EffectiveEIR);
   }
 
-  //Introduce the effect of vaccination. Note that this does not affect cumEIR.
-    expectedNumInfections *= human.getVaccine().getFactor( interventions::Vaccine::PEV );
+  // Only to be consistent with old simulation runs when set to false
+  // Setting this option to true will only affect reporting
+  if(opt_vaccine_genotype == false)
+      //Introduce the effect of vaccination. Note that this does not affect cumEIR.
+      expectedNumInfections *= human.getVaccine().getFactor( interventions::Vaccine::PEV );
 
   //Update pre-erythrocytic immunity
   m_cumulativeEIRa+=effectiveEIR;
@@ -238,8 +243,6 @@ int InfectionIncidenceModel::numNewInfections (Human& human, double effectiveEIR
     if( n > WithinHost::WHInterface::MAX_INFECTIONS ){
         n = WithinHost::WHInterface::MAX_INFECTIONS;
     }
-    mon::reportEventMHI( mon::MHR_NEW_INFECTIONS, human, n );
-    ctsNewInfections += n;
     return n;
   }
   if ( (std::isnan)(expectedNumInfections) ){	// check for not-a-number
@@ -249,4 +252,10 @@ int InfectionIncidenceModel::numNewInfections (Human& human, double effectiveEIR
   return 0;
 }
 
+void InfectionIncidenceModel::reportNumNewInfections(Human& human, int newNumInfections)
+{
+    mon::reportEventMHI( mon::MHR_NEW_INFECTIONS, human, newNumInfections);
+    ctsNewInfections += newNumInfections;
+}
+
 } }

model/Host/InfectionIncidenceModel.h

@@ -107,6 +107,15 @@ class InfectionIncidenceModel
    *        outside transmission due to hospitalisation)
    */
   int numNewInfections(OM::Host::Human& human, double effectiveEIR);
+
+  /** Report the number of new infections
+   * 
+   * This number might be lower than the numNewInfections from this model
+   * due to vaccine factors being applied in the WithinHost model
+   * 
+   * @param int newNumInfections  Number of infections to be reported
+   */
+  void reportNumNewInfections(Human& human, int newNumInfections);
 
 protected:
   /// Calculates the expected number of infections, excluding vaccine effects

model/Host/WithinHost/CommonWithinHost.cpp

@@ -25,6 +25,7 @@
 #include "Host/WithinHost/Genotypes.h"
 #include "Host/WithinHost/Pathogenesis/PathogenesisModel.h"
 #include "util/errors.h"
+#include "util/ModelOptions.h"
 #include "util/AgeGroupInterpolation.h"
 #include "util/random.h"
 #include "util/StreamValidator.h"
@@ -44,7 +45,6 @@ util::AgeGroupInterpolator massByAge;
 
 bool reportInfectedOrPatentInfected = false;
 
-
 // -----  Initialization  -----
 
 void CommonWithinHost::init( const scnXml::Scenario& scenario ){
@@ -69,6 +69,8 @@ CommonWithinHost::CommonWithinHost( LocalRng& rng, double comorbidityFactor ) :
 {
     assert( sim::oneTS() == sim::fromDays(1) || sim::oneTS() == sim::fromDays(5) );
 
+    opt_vaccine_genotype = util::ModelOptions::option (util::VACCINE_GENOTYPE);
+
     // Sample a weight heterogeneity factor
 #ifndef NDEBUG
     int counter = 0;
@@ -135,19 +137,40 @@ void CommonWithinHost::importInfection(LocalRng& rng){
 
 // -----  Density calculations  -----
 
-void CommonWithinHost::update(LocalRng& rng,
-        int nNewInfs, vector<double>& genotype_weights,
-        double ageInYears, double bsvFactor)
+void CommonWithinHost::update(Host::Human &human, LocalRng& rng,
+        int &nNewInfs, vector<double>& genotype_weights,
+        double ageInYears)
 {
     // Note: adding infections at the beginning of the update instead of the end
     // shouldn't be significant since before latentp delay nothing is updated.
-    nNewInfs=min(nNewInfs,MAX_INFECTIONS-numInfs);
-    numInfs += nNewInfs;
+    nNewInfs = min(nNewInfs,MAX_INFECTIONS-numInfs);
+
+    int nNewInfsIgnored = nNewInfs - (MAX_INFECTIONS-numInfs);
+
+    int nNewInfsDiscarded = 0;
+
     assert( numInfs>=0 && numInfs<=MAX_INFECTIONS );
+
     for( int i=0; i<nNewInfs; ++i ) {
         uint32_t genotype = Genotypes::sampleGenotype(rng, genotype_weights);
-        infections.push_back(createInfection (rng, genotype));
+
+        // If opt_vaccine_genotype is true the infection is discarded with probability 1-vaccineFactor
+        if( opt_vaccine_genotype )
+        {
+            double vaccineFactor = human.getVaccine().getFactor( interventions::Vaccine::PEV, genotype );
+            if(vaccineFactor == 1.0 || human.rng().bernoulli(vaccineFactor))
+                infections.push_back(createInfection (rng, genotype));
+            else
+                nNewInfsDiscarded++;
+        }
+        else if (opt_vaccine_genotype == false)
+            infections.push_back(createInfection (rng, genotype));
     }
+
+    // Update nNewInfs, this is the number that will be reported in Human
+    nNewInfs -= nNewInfsDiscarded;
+    numInfs += nNewInfs;
+
     assert( numInfs == static_cast<int>(infections.size()) );
 
     updateImmuneStatus ();
@@ -158,7 +181,6 @@ void CommonWithinHost::update(LocalRng& rng,
 
     bool treatmentLiver = treatExpiryLiver > sim::ts0();
     bool treatmentBlood = treatExpiryBlood > sim::ts0();
-    double survivalFactor_part = bsvFactor * _innateImmSurvFact;
 
     double body_mass = massByAge.eval( ageInYears ) * hetMassMultiplier;
 
@@ -175,7 +197,8 @@ void CommonWithinHost::update(LocalRng& rng,
             if( !expires ){     /* no expiry due to simple treatment model; do update */
                 const double drugFactor = pkpdModel.getDrugFactor(rng, *inf, body_mass);
                 const double immFactor = immunitySurvivalFactor(ageInYears, (*inf)->cumulativeExposureJ());
-                const double survivalFactor = survivalFactor_part * immFactor * drugFactor;
+                const double bsvFactor = human.getVaccine().getFactor(interventions::Vaccine::BSV, opt_vaccine_genotype? (*inf)->genotype() : 0);
+                const double survivalFactor = bsvFactor * _innateImmSurvFact * immFactor * drugFactor;
                 // update, may result in termination of infection:
                 expires = (*inf)->update(rng, survivalFactor, now, body_mass);
             }
@@ -217,6 +240,10 @@ void CommonWithinHost::update(LocalRng& rng,
     for( auto inf = infections.begin(); inf != infections.end(); ++inf ){
         m_y_lag[y_lag_i * Genotypes::N() + (*inf)->genotype()] += (*inf)->getDensity();
     }
+
+    // This is a bug, we keep it this way to be consistent with old simulations
+    if(nNewInfsIgnored > 0)
+        nNewInfs += nNewInfsIgnored;
 }
 
 void CommonWithinHost::addProphylacticEffects(const vector<double>& pClearanceByTime) {

model/Host/WithinHost/CommonWithinHost.h

@@ -51,8 +51,8 @@ class CommonWithinHost : public WHFalciparum
     virtual void treatPkPd(size_t schedule, size_t dosage, double age, double delay_d);
     virtual void clearImmunity();
 
-    virtual void update (LocalRng& rng, int nNewInfs, vector<double>& genotype_weights,
-            double ageInYears, double bsvFactor);
+    virtual void update (Host::Human &human, LocalRng& rng, int &nNewInfs, vector<double>& genotype_weights,
+            double ageInYears);
 
     virtual void addProphylacticEffects(const vector<double>& pClearanceByTime);
 
@@ -88,6 +88,8 @@ class CommonWithinHost : public WHFalciparum
      * the idea is that each WithinHostModel has its own list of infections. */
     //TODO: better to template class over infection type than use dynamic type?
     std::list<CommonInfection*> infections;
+
+    bool opt_vaccine_genotype = false;
 };
 
 } }

model/Host/WithinHost/DescriptiveWithinHost.cpp

@@ -37,7 +37,6 @@ namespace WithinHost {
 
 extern bool bugfix_max_dens;    // DescriptiveInfection.cpp
 bool reportPatentInfected = false;
-
 // -----  Initialization  -----
 
 void DescriptiveWithinHostModel::initDescriptive(){
@@ -48,6 +47,7 @@ DescriptiveWithinHostModel::DescriptiveWithinHostModel( LocalRng& rng, double co
         WHFalciparum( rng, comorbidityFactor )
 {
     assert( sim::oneTS() == sim::fromDays(5) );
+    opt_vaccine_genotype = util::ModelOptions::option (util::VACCINE_GENOTYPE);
 }
 
 DescriptiveWithinHostModel::~DescriptiveWithinHostModel() {}
@@ -98,18 +98,30 @@ void DescriptiveWithinHostModel::importInfection(LocalRng& rng){
 
 // -----  Density calculations  -----
 
-void DescriptiveWithinHostModel::update(LocalRng& rng,
-        int nNewInfs, vector<double>& genotype_weights,
-        double ageInYears, double bsvFactor)
+void DescriptiveWithinHostModel::update(Host::Human &human, LocalRng& rng,
+        int &nNewInfs, vector<double>& genotype_weights,
+        double ageInYears)
 {
     // Note: adding infections at the beginning of the update instead of the end
     // shouldn't be significant since before latentp delay nothing is updated.
-    nNewInfs=min(nNewInfs,MAX_INFECTIONS-numInfs);
+    int nNewInfsToBeCreated = nNewInfs;
+
+    nNewInfs = min(nNewInfs,MAX_INFECTIONS-numInfs);
+
     numInfs += nNewInfs;
     assert( numInfs>=0 && numInfs<=MAX_INFECTIONS );
     for( int i=0; i<nNewInfs; ++i ) {
         uint32_t genotype = Genotypes::sampleGenotype(rng, genotype_weights);
-        infections.push_back(DescriptiveInfection (rng, genotype));
+
+        // If opt_vaccine_genotype is true the infection is discarded with probability 1-vaccineFactor
+        if( opt_vaccine_genotype )
+        {
+            double vaccineFactor = human.getVaccine().getFactor( interventions::Vaccine::PEV, genotype );
+            if(vaccineFactor == 1.0 || human.rng().bernoulli(vaccineFactor))
+                infections.push_back(DescriptiveInfection (rng, genotype));
+        }
+        else if (opt_vaccine_genotype == false)
+            infections.push_back(DescriptiveInfection (rng, genotype));
     }
     assert( numInfs == static_cast<int>(infections.size()) );
 
@@ -142,6 +154,8 @@ void DescriptiveWithinHostModel::update(LocalRng& rng,
         // See MAX_DENS_CORRECTION in DescriptiveInfection.cpp.
         double infStepMaxDens = timeStepMaxDensity;
         double immSurvFact = immunitySurvivalFactor(ageInYears, inf->cumulativeExposureJ());
+        double bsvFactor = human.getVaccine().getFactor(interventions::Vaccine::BSV, opt_vaccine_genotype? inf->genotype() : 0);
+
         inf->determineDensities(rng, m_cumulative_h, infStepMaxDens, immSurvFact, _innateImmSurvFact, bsvFactor);
 
         if (bugfix_max_dens)
@@ -172,6 +186,10 @@ void DescriptiveWithinHostModel::update(LocalRng& rng,
     for( auto inf = infections.begin(); inf != infections.end(); ++inf ){
         m_y_lag[y_lag_i * Genotypes::N() + inf->genotype()] += inf->getDensity();
     }
+
+    // This is a bug, we keep it this way to be consistent with old simulations
+    if(opt_vaccine_genotype == false)
+        nNewInfs = nNewInfsToBeCreated;
 }
 
 

model/Host/WithinHost/DescriptiveWithinHost.h

@@ -49,8 +49,8 @@ class DescriptiveWithinHostModel : public WHFalciparum {
     virtual void loadInfection(istream& stream);
     virtual void clearImmunity();
 
-    virtual void update(LocalRng& rng, int nNewInfs, vector<double>& genotype_weights,
-            double ageInYears, double bsvFactor);
+    virtual void update(Host::Human &human, LocalRng& rng, int &nNewInfs, vector<double>& genotype_weights,
+            double ageInYears);
 
     virtual bool summarize( Host::Human& human )const;
 
@@ -68,6 +68,8 @@ class DescriptiveWithinHostModel : public WHFalciparum {
      * Since infection models and within host models are very much intertwined,
      * the idea is that each WithinHostModel has its own list of infections. */
      std::list<DescriptiveInfection> infections;
+
+     bool opt_vaccine_genotype = false;
 };
 
 } }

model/Host/WithinHost/WHFalciparum.cpp

@@ -172,7 +172,7 @@ const double PTM_tau= 0.066;
 const double PTM_tau_prime = 1.0 / sqrt(1.0 / PTM_tau);
 const double PTM_mu= -8.1;
 
-double WHFalciparum::probTransmissionToMosquito( double tbvFactor, double *sumX ) const{
+double WHFalciparum::probTransmissionToMosquito( double *sumX ) const{
     // This model (often referred to as the gametocyte model) was designed for
     // 5-day time steps. We use the same model (sampling 10, 15 and 20 days
     // ago) for 1-day time steps to avoid having to design and analyse a new
@@ -206,7 +206,6 @@ double WHFalciparum::probTransmissionToMosquito( double tbvFactor, double *sumX
     pTransmit=std::min(pTransmit, 1.0);
 
     // Include here the effect of transmission-blocking vaccination:
-    pTransmit *= tbvFactor;
     util::streamValidate( pTransmit );
     return pTransmit;
 }

model/Host/WithinHost/WHFalciparum.h

@@ -57,7 +57,7 @@ class WHFalciparum : public WHInterface {
     virtual ~WHFalciparum();
     //@}
 
-    virtual double probTransmissionToMosquito( double tbvFactor, double *sumX )const;
+    virtual double probTransmissionToMosquito( double *sumX )const;
     virtual double pTransGenotype( double pTrans, double sumX, size_t genotype );
 
     // No PQ treatment for falciparum in current models:

model/Host/WithinHost/WHInterface.h

@@ -113,8 +113,8 @@ class WHInterface {
      * 
      * Calculates the value during the call, which is expensive (cache externally
      * if the value is needed multiple times). */
-    virtual double probTransmissionToMosquito( double tbvFactor,
-                                               double *sumX )const =0;
+    virtual double probTransmissionToMosquito(double *sumX )const =0;
+
     /** Calculates a probability of transmitting an infection of a given
      * genotype to a mosquito, given the two outputs of
      * probTransmissionToMosquito(). Only available for WHFalciparum and
@@ -166,8 +166,8 @@ class WHInterface {
      * @param ageInYears Age of human
      * @param bsvFactor Parasite survival factor for blood-stage vaccines
      */
-    virtual void update(LocalRng& rng, int nNewInfs, vector<double>& genotype_weights,
-            double ageInYears, double bsvFactor) =0;
+    virtual void update(Host::Human &human, LocalRng& rng, int &nNewInfs, vector<double>& genotype_weights,
+            double ageInYears) =0;
 
     /** TODO: this should not need to be exposed. It is currently used by a
      * severe outcome (pDeath) model inside the EventScheduler "case

model/Host/WithinHost/WHVivax.cpp

@@ -311,14 +311,14 @@ WHVivax::~WHVivax(){
 #endif
 }
 
-double WHVivax::probTransmissionToMosquito( double tbvFactor, double *sumX )const{
+double WHVivax::probTransmissionToMosquito( double *sumX )const{
     assert( WithinHost::Genotypes::N() == 1 );
     for(auto inf = infections.begin();
          inf != infections.end(); ++inf)
     {
         if( inf->isPatent() ){
             // we have gametocytes from at least one brood
-            return probBloodStageInfectiousToMosq * tbvFactor;
+            return probBloodStageInfectiousToMosq;
         }
     }
     return 0;   // no gametocytes
@@ -349,9 +349,9 @@ void WHVivax::importInfection(LocalRng& rng){
     infections.push_back( VivaxBrood( rng, this ) );
 }
 
-void WHVivax::update(LocalRng& rng,
-        int nNewInfs, vector<double>&,
-        double ageInYears, double)
+void WHVivax::update(Host::Human &human, LocalRng& rng,
+        int &nNewInfs, vector<double>&,
+        double ageInYears)
 {
     pSevere = 0.0;
 

model/Host/WithinHost/WHVivax.h

@@ -147,15 +147,15 @@ class WHVivax : public WHInterface {
     virtual ~WHVivax();
     //@}
 
-    virtual double probTransmissionToMosquito( double tbvFactor, double *sumX )const;
+    virtual double probTransmissionToMosquito( double *sumX )const;
     virtual double pTransGenotype( double pTrans, double sumX, size_t genotype );
 
     virtual bool summarize(Host::Human& human) const;
 
     virtual void importInfection(LocalRng& rng);
 
-    virtual void update(LocalRng& rng, int nNewInfs, vector<double>& genotype_weights,
-            double ageInYears, double bsvFactor);
+    virtual void update(Host::Human &human, LocalRng& rng, int &nNewInfs, vector<double>& genotype_weights,
+            double ageInYears);
 
     virtual bool diagnosticResult( LocalRng& rng, const Diagnostic& diagnostic ) const;