Merge pull request #1058 from amath-idm/rc3.0.4

Rc3.0.4

CHANGELOG.rst

@@ -24,6 +24,18 @@ These are the major improvements we are currently working on. If there is a spec
 Latest versions (3.0.x)
 ~~~~~~~~~~~~~~~~~~~~~~~
 
+Version 3.0.4 (2021-05-19)
+--------------------------
+- Fixed a bug that prevented simulations from being run *without* prognoses by age.
+- Fixed an array length mismatch for single-dose vaccines.
+- The default antibody kinetics are now a 3-part curve, with a 14-day growth, 250 day exp decay and then another exponential decay with a exponentially decaying decay parameter. This is captured in the new NAb functional form, ``nab_growth_decay``. To align with this change, NAbs are now initialized at the time of infection, so that individuals build immunity over the course of infection. 
+- Some strain parameter changes based on https://www.eurosurveillance.org/content/10.2807/1560-7917.ES.2021.26.16.2100348
+- Added strain to the infection log
+- Removed the ``rel_imm_strain`` parameter; self-immunity is now always 1.0.
+- Updated vaccine and strain parameter values based on fits to empirical data.
+- Merged multisims now use the labels from each multisim, rather than the sim labels, for plotting.
+- *GitHub info*: PR `1058 <https://github.com/amath-idm/covasim/pull/1058>`__
+
 
 Version 3.0.3 (2021-05-17)
 --------------------------

covasim/analysis.py

@@ -1266,6 +1266,8 @@ class Calibration(Analyzer):
         calib = cv.Calibration(sim, calib_pars, total_trials=100)
         calib.calibrate()
         calib.plot()
+
+    New in version 3.0.3.
     '''
 
     def __init__(self, sim, calib_pars=None, custom_fn=None, n_trials=None, n_workers=None, total_trials=None, name=None, db_name=None, storage=None, label=None, verbose=True):

covasim/defaults.py

@@ -94,7 +94,8 @@ def __init__(self):
         # Neutralizing antibody states, not by strain
         self.nab_states = [
             'prior_symptoms',   # Float
-            'init_nab',         # Float, initial neutralization titre relative to convalescent plasma
+            'peak_nab',         # Float, peak neutralization titre relative to convalescent plasma
+            'last_nab',         # Float, neutralization titre relative to convalescent plasma on day of boost
             'nab',              # Float, current neutralization titre relative to convalescent plasma
         ]
 

covasim/immunity.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import sciris as sc
+from scipy.special import expit
 from . import utils as cvu
 from . import defaults as cvd
 from . import parameters as cvpar
@@ -146,7 +147,7 @@ def get_vaccine_pars(pars):
 
 def init_nab(people, inds, prior_inf=True):
     '''
-    Draws an initial neutralizing antibody (NAb) level for individuals.
+    Draws a peak neutralizing antibody (NAb) level for individuals.
     Can come from a natural infection or vaccination and depends on if there is prior immunity:
     1) a natural infection. If individual has no existing NAb, draw from distribution
     depending upon symptoms. If individual has existing NAb, multiply booster impact
@@ -157,7 +158,7 @@ def init_nab(people, inds, prior_inf=True):
     nab_arrays = people.nab[inds]
     prior_nab_inds = cvu.idefined(nab_arrays, inds) # Find people with prior NAb
     no_prior_nab_inds = np.setdiff1d(inds, prior_nab_inds) # Find people without prior NAb
-    peak_nab = people.init_nab[prior_nab_inds]
+    peak_nab = people.peak_nab[prior_nab_inds] # Find the prior peak for those with prior NAbs
     pars = people.pars
 
     # NAb from infection
@@ -168,12 +169,14 @@ def init_nab(people, inds, prior_inf=True):
             init_nab = cvu.sample(**pars['nab_init'], size=len(no_prior_nab_inds))
             prior_symp = people.prior_symptoms[no_prior_nab_inds]
             no_prior_nab = (2**init_nab) * prior_symp
-            people.init_nab[no_prior_nab_inds] = no_prior_nab
+            people.peak_nab[no_prior_nab_inds] = no_prior_nab
 
         # 2) Prior NAb: multiply existing NAb by boost factor
         if len(prior_nab_inds):
+            last_nab = people.nab[prior_nab_inds]
+            people.last_nab[prior_nab_inds] = last_nab
             init_nab = peak_nab * nab_boost
-            people.init_nab[prior_nab_inds] = init_nab
+            people.peak_nab[prior_nab_inds] = init_nab
 
     # NAb from a vaccine
     else:
@@ -182,33 +185,45 @@ def init_nab(people, inds, prior_inf=True):
         # 1) No prior NAb: draw NAb from a distribution and compute
         if len(no_prior_nab_inds):
             init_nab = cvu.sample(**vaccine_pars['nab_init'], size=len(no_prior_nab_inds))
-            people.init_nab[no_prior_nab_inds] = 2**init_nab
+            people.peak_nab[no_prior_nab_inds] = 2**init_nab
 
         # 2) Prior nab (from natural or vaccine dose 1): multiply existing nab by boost factor
         if len(prior_nab_inds):
+            last_nab = people.nab[prior_nab_inds]
+            people.last_nab[prior_nab_inds] = last_nab
             nab_boost = vaccine_pars['nab_boost']  # Boosting factor for vaccination
             init_nab = peak_nab * nab_boost
-            people.init_nab[prior_nab_inds] = init_nab
+            people.peak_nab[prior_nab_inds] = init_nab
 
     return
 
 
 def check_nab(t, people, inds=None):
-    ''' Determines current NAb based on date since recovered/vaccinated.'''
-
+    ''' Determines current NAb based on date since recovered/vaccinated.
+        First step: determine if we are in the growth or decay period
+            If in growth, pull nabs of inds and add % of peak nabs
+            If in decay, % of peak nabs
+    '''
     # Indices of people who've had some nab event
-    rec_inds = cvu.defined(people.date_recovered[inds])
+    inf_inds = cvu.defined(people.date_exposed[inds])
     vac_inds = cvu.defined(people.date_vaccinated[inds])
-    both_inds = np.intersect1d(rec_inds, vac_inds)
+    both_inds = np.intersect1d(inf_inds, vac_inds)
 
     # Time since boost
     t_since_boost = np.full(len(inds), np.nan, dtype=cvd.default_int)
-    t_since_boost[rec_inds] = t-people.date_recovered[inds[rec_inds]]
+    t_since_boost[inf_inds] = t-people.date_exposed[inds[inf_inds]]
     t_since_boost[vac_inds] = t-people.date_vaccinated[inds[vac_inds]]
-    t_since_boost[both_inds] = t-np.maximum(people.date_recovered[inds[both_inds]],people.date_vaccinated[inds[both_inds]])
+    t_since_boost[both_inds] = t-np.maximum(people.date_exposed[inds[both_inds]],people.date_vaccinated[inds[both_inds]])
+
+    # Determine which nabs are in decay (peak > current)
+    nabs = people.last_nab[inds]
+    if people.pars['nab_decay']['form'] == 'nab_growth_decay':
+        inds_in_decay = cvu.true(t_since_boost >= people.pars['nab_decay']['growth_time'])
+        nabs[inds_in_decay] = 0
+    nabs = np.nan_to_num(nabs)
 
     # Set current NAb
-    people.nab[inds] = people.pars['nab_kin'][t_since_boost] * people.init_nab[inds]
+    people.nab[inds] = nabs + people.pars['nab_kin'][t_since_boost] * people.peak_nab[inds]
 
     return
 
@@ -234,7 +249,9 @@ def nab_to_efficacy(nab, ax, function_args):
     if ax == 'sus':
         slope = args['slope']
         n_50 = args['n_50']
-        efficacy = 1 / (1 + np.exp(-slope * (np.log10(nab) - np.log10(n_50))))  # from logistic regression computed in R using data from Khoury et al
+        logNAb = np.log(nab)
+        VE = logNAb * slope + n_50
+        efficacy = expit(VE)  # from linear in logit-log space (see fit https://github.com/amath-idm/COVID-Immune-Modeling/blob/main/scripts/nab2eff.R)
     else:
         efficacy = np.full(len(nab), fill_value=args)
     return efficacy
@@ -251,7 +268,7 @@ def init_immunity(sim, create=False):
         return
 
     # Pull out all of the circulating strains for cross-immunity
-    ns       = sim['n_strains']
+    ns = sim['n_strains']
 
     # If immunity values have been provided, process them
     if sim['immunity'] is None or create:
@@ -375,14 +392,18 @@ def precompute_waning(length, pars=None):
     pars = sc.dcp(pars)
     form = pars.pop('form')
     choices = [
-        'nab_decay', # Default if no form is provided
+        'nab_growth_decay', # Default if no form is provided
+        'nab_decay',
         'exp_decay',
         'linear_growth',
         'linear_decay'
     ]
 
     # Process inputs
-    if form is None or form == 'nab_decay':
+    if form is None or form == 'nab_growth_decay':
+        output = nab_growth_decay(length, **pars)
+
+    elif form == 'nab_decay':
         output = nab_decay(length, **pars)
 
     elif form == 'exp_decay':
@@ -402,6 +423,40 @@ def precompute_waning(length, pars=None):
     return output
 
 
+def nab_growth_decay(length, growth_time, decay_rate1, decay_time1, decay_rate2):
+    '''
+    Returns an array of length 'length' containing the evaluated function nab growth/decay
+    function at each point.
+    Uses linear growth + exponential decay, with the rate of exponential decay also set to
+    exponentially decay (!) after 250 days.
+    Args:
+        length (int): number of points
+        growth_time (int): length of time NAbs grow (used to determine slope)
+        decay_rate1 (float): initial rate of exponential decay
+        decay_time1 (float): time on the first exponential decay
+        decay_rate2 (float): the rate at which the decay decays
+    '''
+
+    def f1(t, growth_time):
+        '''Simple linear growth'''
+        return (1 / growth_time) * t
+
+    def f2(t, decay_rate1):
+        ''' Simple exponential decay '''
+        return np.exp(-t*decay_rate1)
+
+    def f3(t, decay_rate1, decay_time1, decay_rate2):
+        ''' Complex exponential decay '''
+        return np.exp(-t*(decay_rate1*np.exp(-(t-decay_time1)*decay_rate2)))
+
+    t  = np.arange(length, dtype=cvd.default_int)
+    y1 = f1(cvu.true(t<= growth_time), growth_time)
+    y2 = f2(cvu.true((t>growth_time)*(t<= decay_time1)), decay_rate1)
+    y3 = f3(cvu.true(t>decay_time1), decay_rate1, decay_time1, decay_rate2)
+    y  = np.concatenate([y1, y2, y3])
+
+    return y
+
 def nab_decay(length, decay_rate1, decay_time1, decay_rate2):
     '''
     Returns an array of length 'length' containing the evaluated function nab decay

covasim/interventions.py

@@ -1324,8 +1324,6 @@ def initialize(self, sim):
 
 
         self.days = process_days(sim, self.days) # days that group becomes eligible
-        self.first_dose_nab_days  = [None]*sim.npts # People who get nabs from first dose
-        self.second_dose_nab_days = [None]*sim.npts # People who get nabs from second dose (if relevant)
         self.second_dose_days     = [None]*sim.npts # People who get second dose (if relevant)
         self.vaccinated           = [None]*sim.npts # Keep track of inds of people vaccinated on each day
         self.vaccinations         = np.zeros(sim['pop_size'], dtype=cvd.default_int) # Number of doses given per person
@@ -1344,7 +1342,7 @@ def apply(self, sim):
 
             # Vaccinate people with their first dose
             vacc_inds = np.array([], dtype=int) # Initialize in case no one gets their first dose
-            for ind in find_day(self.days, sim.t, interv=self, sim=sim):
+            for _ in find_day(self.days, sim.t, interv=self, sim=sim):
                 vacc_probs = np.zeros(sim['pop_size'])
                 unvacc_inds = sc.findinds(~sim.people.vaccinated)
                 if self.subtarget is not None:
@@ -1364,16 +1362,10 @@ def apply(self, sim):
                         next_dose_day = sim.t + self.p.interval
                         if next_dose_day < sim['n_days']:
                             self.second_dose_days[next_dose_day] = vacc_inds
-                            self.first_dose_nab_days[next_dose_day] = vacc_inds
-                    else:
-                        self.first_dose_nab_days[sim.t] = vacc_inds
 
             # Also, if appropriate, vaccinate people with their second dose
             vacc_inds_dose2 = self.second_dose_days[sim.t]
             if vacc_inds_dose2 is not None:
-                next_nab_day = sim.t + self.p.interval
-                if next_nab_day < sim['n_days']:
-                    self.second_dose_nab_days[next_nab_day] = vacc_inds_dose2
                 vacc_inds = np.concatenate((vacc_inds, vacc_inds_dose2), axis=None)
                 sim.people.flows['new_vaccinations'] += len(vacc_inds_dose2)
 
@@ -1386,17 +1378,8 @@ def apply(self, sim):
 
                 # Update vaccine attributes in sim
                 sim.people.vaccinations[vacc_inds] = self.vaccinations[vacc_inds]
-
-            # check whose nabs kick in today and then init_nabs for them!
-            vaccine_nabs_first_dose_inds = self.first_dose_nab_days[sim.t]
-            vaccine_nabs_second_dose_inds = self.second_dose_nab_days[sim.t]
-
-            vaccine_nabs_inds = [vaccine_nabs_first_dose_inds, vaccine_nabs_second_dose_inds]
-
-            for vaccinees in vaccine_nabs_inds:
-                if vaccinees is not None:
-                    sim.people.date_vaccinated[vaccinees] = sim.t
-                    cvi.init_nab(sim.people, vaccinees, prior_inf=False)
+                sim.people.date_vaccinated[vacc_inds] = sim.t
+                cvi.init_nab(sim.people, vacc_inds, prior_inf=False)
 
         return
 

covasim/parameters.py

@@ -68,10 +68,10 @@ def make_pars(set_prognoses=False, prog_by_age=True, version=None, **kwargs):
     # Parameters used to calculate immunity
     pars['use_waning']      = False # Whether to use dynamically calculated immunity
     pars['nab_init']        = dict(dist='normal', par1=0, par2=2)  # Parameters for the distribution of the initial level of log2(nab) following natural infection, taken from fig1b of https://doi.org/10.1101/2021.03.09.21252641
-    pars['nab_decay']       = dict(form='nab_decay', decay_rate1=np.log(2)/90, decay_time1=250, decay_rate2=0.001) # Parameters describing the kinetics of decay of nabs over time, taken from fig3b of https://doi.org/10.1101/2021.03.09.21252641
+    pars['nab_decay']       = dict(form='nab_growth_decay', growth_time=14, decay_rate1=np.log(2) / 90, decay_time1=250, decay_rate2=0.001)
     pars['nab_kin']         = None # Constructed during sim initialization using the nab_decay parameters
     pars['nab_boost']       = 1.5 # Multiplicative factor applied to a person's nab levels if they get reinfected. # TODO: add source
-    pars['nab_eff']         = dict(sus=dict(slope=1.6, n_50=0.05), symp=0.1, sev=0.52) # Parameters to map nabs to efficacy
+    pars['nab_eff']         = dict(sus=dict(slope=0.8412132, n_50=1.318724), symp=0.3, sev=0.52) # Parameters to map nabs to efficacy
     pars['rel_imm_symp']    = dict(asymp=0.85, mild=1, severe=1.5) # Relative immunity from natural infection varies by symptoms
     pars['immunity']        = None  # Matrix of immunity and cross-immunity factors, set by init_immunity() in immunity.py
 
@@ -367,23 +367,23 @@ def get_strain_pars(default=False):
             rel_symp_prob   = 1.0, # Inconclusive evidence on the likelihood of symptom development. See https://www.thelancet.com/journals/lanpub/article/PIIS2468-2667(21)00055-4/fulltext
             rel_severe_prob = 1.8, # From https://www.ssi.dk/aktuelt/nyheder/2021/b117-kan-fore-til-flere-indlaggelser and https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/961042/S1095_NERVTAG_update_note_on_B.1.1.7_severity_20210211.pdf
             rel_crit_prob   = 1.0, # Various studies have found increased mortality for B117 (summary here: https://www.thelancet.com/journals/laninf/article/PIIS1473-3099(21)00201-2/fulltext#tbl1), but not necessarily when conditioned on having developed severe disease
-            rel_death_prob  = 1.0, # See comment above.
+            rel_death_prob  = 1.6, # See comment above.
         ),
 
         b1351 = dict(
             rel_beta        = 1.4,
-            rel_symp_prob   = 1.0,
-            rel_severe_prob = 1.4,
+            rel_symp_prob   = 2.0,
+            rel_severe_prob = 3.6,
             rel_crit_prob   = 1.0,
-            rel_death_prob  = 1.4,
+            rel_death_prob  = 3.2,
         ),
 
         p1 = dict(
             rel_beta        = 1.4, # Estimated to be 1.7–2.4-fold more transmissible than wild-type: https://science.sciencemag.org/content/early/2021/04/13/science.abh2644
             rel_symp_prob   = 1.0,
-            rel_severe_prob = 1.4,
+            rel_severe_prob = 2.6,
             rel_crit_prob   = 1.0,
-            rel_death_prob  = 2.0,
+            rel_death_prob  = 1.67,
         )
     )
 
@@ -496,47 +496,47 @@ def get_vaccine_dose_pars(default=False):
     pars = dict(
 
         default = dict(
-            nab_eff   = dict(sus=dict(slope=1.6, n_50=0.05)),
+            nab_eff   = dict(sus=dict(slope=0.8412132, n_50=1.318724)),
             nab_init  = dict(dist='normal', par1=2, par2=2),
             nab_boost = 2,
             doses     = 1,
             interval  = None,
         ),
 
         pfizer = dict(
-            nab_eff   = dict(sus=dict(slope=1.6, n_50=0.05)),
+            nab_eff   = dict(sus=dict(slope=0.8412132, n_50=1.318724)),
             nab_init  = dict(dist='normal', par1=2, par2=2),
             nab_boost = 3,
             doses     = 2,
             interval  = 21,
         ),
 
         moderna = dict(
-            nab_eff   = dict(sus=dict(slope=1.6, n_50=0.05)),
+            nab_eff   = dict(sus=dict(slope=0.8412132, n_50=1.318724)),
             nab_init  = dict(dist='normal', par1=2, par2=2),
             nab_boost = 3,
             doses     = 2,
             interval  = 28,
         ),
 
         az = dict(
-            nab_eff   = dict(sus=dict(slope=1.6, n_50=0.05)),
-            nab_init  = dict(dist='normal', par1=-0.85, par2=2),
+            nab_eff   = dict(sus=dict(slope=0.8412132, n_50=1.318724)),
+            nab_init  = dict(dist='normal', par1=-1, par2=2),
             nab_boost = 3,
             doses     = 2,
             interval  = 21,
         ),
 
         jj = dict(
-            nab_eff   = dict(sus=dict(slope=1.6, n_50=0.05)),
-            nab_init  = dict(dist='normal', par1=-1.1, par2=2),
+            nab_eff   = dict(sus=dict(slope=0.8412132, n_50=1.318724)),
+            nab_init  = dict(dist='normal', par1=1, par2=2),
             nab_boost = 3,
             doses     = 1,
             interval  = None,
         ),
 
         novavax = dict(
-            nab_eff   = dict(sus=dict(slope=1.6, n_50=0.05)),
+            nab_eff   = dict(sus=dict(slope=0.8412132, n_50=1.318724)),
             nab_init  = dict(dist='normal', par1=-0.9, par2=2),
             nab_boost = 3,
             doses     = 2,