Merge pull request #260 from SABS-R3-Epidemiology/waning-immunity

Rate multipliers and IgG count

README.md

@@ -86,6 +86,7 @@ Configure a simulation with a number of parameters. These are split into three c
 * `initial_infected_number`: The initial number of infected individuals in the population
 * `initial_infect_cell`: Whether to choose initial infected individuals from a single cell
 * `simulation_seed`:  Random seed for reproducible simulations - see above _(Optional)_
+* `include_waning`: Boolean to determine whether immunity waning is included in the simulation _(Default false)_
 
 *`file_params`* _(For controlling output location)_
 * `output_file`: String for the name of the output .csv file

pyEpiabm/README.rst

@@ -108,6 +108,7 @@ sim_params
 - ``initial_infected_number``: The initial number of infected individuals in the population
 - ``initial_infect_cell``: Whether to choose initial infected individuals from a single cell
 - ``simulation_seed``: Random seed for reproducible simulations - see above *(Optional)*
+- ``include_waning``: Boolean to determine whether immunity waning is included in the simulation *(Default false)*
 
 file_params
 """""""""""
@@ -123,9 +124,9 @@ inf_history_params
 *(For controlling the infection history output - Default None)*
 
 - ``output_dir``: String for the location for the output files, as a relative path
-- ``status_output``: Boolean to determine whether we need a csv file containing infection status values (Default false)
-- ``infectiousness_output``: Boolean to determine whether we need a csv file containing infectiousness (viral load) values (Default false)
-- ``compress``: Boolean to determine whether we compress a csv file containing infection status values and/or a csv file containing infectiousness (viral load) values if they are written (Default false)
+- ``status_output``: Boolean to determine whether we need a csv file containing infection status values *(Default false)*
+- ``infectiousness_output``: Boolean to determine whether we need a csv file containing infectiousness (viral load) values *(Default false)*
+- ``compress``: Boolean to determine whether we compress a csv file containing infection status values and/or a csv file containing infectiousness (viral load) values if they are written *(Default false)*
 
 Two lists of sweeps must also be passed to this function - the first
 will be executed once at the start of the simulation (i.e. to determine

pyEpiabm/docs/source/utility.rst

@@ -8,11 +8,17 @@ Utility provides various methods that don't act on a population directly.
 
 Overview:
 
+- :class:`AntibodyMultiplier`
 - :class:`DistanceFunctions`
 - :class:`InverseCdf`
 - :class:`RandomMethods`
+- :class:`RateMultiplier`
 - :class:`SpatialKernel`
 
+.. autoclass:: AntibodyMultiplier
+    :members:
+    :special-members: __call__
+
 .. autoclass:: DistanceFunctions
     :members:
 
@@ -22,6 +28,10 @@ Overview:
 .. autoclass:: RandomMethods
     :members:
 
+.. autoclass:: RateMultiplier
+    :members:
+    :special-members: __call__
+
 .. autoclass:: SpatialKernel
     :members:
 

pyEpiabm/pyEpiabm/core/person.py

@@ -52,6 +52,7 @@ def __init__(self, microcell, age_group=None):
         self.time_of_status_change = None
         self.infection_start_time = None
         self.time_of_recovery = None
+        self.num_times_infected = 0
         self.care_home_resident = False
         self.key_worker = False
         self.date_positive = None
@@ -146,9 +147,10 @@ def update_status(self,
         if self.infection_status == InfectionStatus.Susceptible and \
                 self.household is not None:
             self.household.add_susceptible_person(self)
-        if self.infection_status == InfectionStatus.Exposed and \
-                self.household is not None:
-            self.household.remove_susceptible_person(self)
+        if self.infection_status == InfectionStatus.Exposed:
+            self.increment_num_times_infected()
+            if self.household is not None:
+                self.household.remove_susceptible_person(self)
 
     def add_place(self, place, person_group: int = 0):
         """Method adds a place to the place list if the person visits
@@ -270,3 +272,9 @@ def set_time_of_recovery(self, time: float):
 
         """
         self.time_of_recovery = time
+
+    def increment_num_times_infected(self):
+        """Increments the number of times the person has been Infected as a
+        useful parameter to keep track of.
+        """
+        self.num_times_infected += 1

pyEpiabm/pyEpiabm/output/_csv_dict_writer.py

@@ -68,4 +68,5 @@ def compress(self):
         logging.info(f"Zip file created for {self.filename}")
         df = pd.read_csv(self.filepath)
         df.to_csv(output_filepath, index=False, compression={'method': 'zip'})
+        self.f.close()
         os.remove(self.filepath)

pyEpiabm/pyEpiabm/property/household_foi.py

@@ -65,8 +65,8 @@ def household_susc(infector, infectee, time: float):
             Susceptibility parameter of household
 
         """
-        household_susceptibility = PersonalInfection.person_susc(
-            infector, infectee, time)
+        household_susceptibility = PersonalInfection.person_susc(infectee,
+                                                                 time)
         if (hasattr(infectee.microcell, 'distancing_start_time')) and (
                 infectee.microcell.distancing_start_time is not None) and (
                     infectee.microcell.distancing_start_time <= time):

pyEpiabm/pyEpiabm/property/personal_foi.py

@@ -1,8 +1,10 @@
 #
 # Calculate infectiousness and susceptibility for an individual
 #
+from collections import defaultdict
 
 from pyEpiabm.core import Parameters
+from pyEpiabm.utility import AntibodyMultiplier
 
 
 class PersonalInfection:
@@ -37,16 +39,14 @@ def person_inf(infector, time: float):
         return infector_inf
 
     @staticmethod
-    def person_susc(infector, infectee, time: float):
+    def person_susc(infectee, time: float):
         """Calculate the susceptibility of one person to another.
 
         Does not yet import WAIFW matrix from Polymod data to determine
         age dependant contact between individuals.
 
         Parameters
         ----------
-        infector : Person
-            Infector
         infectee : Person
             Infectee
         time : float
@@ -55,8 +55,24 @@ def person_susc(infector, infectee, time: float):
         Returns
         -------
         float
-            Susceptibility parameter of household
+            Susceptibility parameter of person
 
         """
-
-        return 1.0
+        # If we are using waning immunity then we use a multiplier from
+        # igg_foi_multiplier. Otherwise, we set the susceptibility to 1.0.
+        if Parameters.instance().use_waning_immunity and\
+                (infectee.num_times_infected >= 1):
+            params = defaultdict(int,
+                                 Parameters.instance().antibody_level_params)
+            if not hasattr(PersonalInfection, 'm'):
+                PersonalInfection.m =\
+                    AntibodyMultiplier(params['igg_peak_at_age_41'],
+                                       params['igg_half_life_at_age_41'],
+                                       params['peak_change_per_10_yrs_age'],
+                                       params['half_life_diff_per_10_yrs_age'],
+                                       params['days_positive_pcr_to_max_igg'])
+            time_since_infection = time - infectee.infection_start_time
+            return 1.0 * PersonalInfection.m(time_since_infection,
+                                             infectee.age_group)
+        else:
+            return 1.0

pyEpiabm/pyEpiabm/property/place_foi.py

@@ -76,7 +76,7 @@ def place_susc(place, infectee, time: float):
             Susceptibility parameter of place
 
         """
-        place_susc = 1.0
+        place_susc = PersonalInfection.person_susc(infectee, time)
         place_idx = place.place_type.value - 1
         if (hasattr(infectee.microcell, 'distancing_start_time')) and (
                 infectee.microcell.distancing_start_time is not None) and (

pyEpiabm/pyEpiabm/property/spatial_foi.py

@@ -7,6 +7,8 @@
 import pyEpiabm.core
 from pyEpiabm.core import Parameters
 
+from .personal_foi import PersonalInfection
+
 
 class SpatialInfection:
     """Class to calculate the infectiousness and susceptibility
@@ -106,9 +108,9 @@ def spatial_susc(susc_cell, infectee, time: float):
             Susceptibility parameter of cell
 
         """
-        spatial_susc = 1.0
+        spatial_susc = PersonalInfection.person_susc(infectee, time)
         if pyEpiabm.core.Parameters.instance().use_ages:
-            spatial_susc = pyEpiabm.core.Parameters.instance().\
+            spatial_susc *= pyEpiabm.core.Parameters.instance().\
                 age_contact[infectee.age_group]
 
         spatial_susc *= Parameters.instance().\

pyEpiabm/pyEpiabm/routine/file_population_config.py

@@ -125,6 +125,7 @@ def make_pop(input_file: str, random_seed: int = None, time: float = 0):
                         if str(person.infection_status).startswith('Infect'):
                             HostProgressionSweep.set_infectiousness(person,
                                                                     time)
+                            person.increment_num_times_infected()
 
             # Add households and places to microcell
             if len(Parameters.instance().household_size_distribution) == 0:

pyEpiabm/pyEpiabm/sweep/host_progression_sweep.py

@@ -2,6 +2,8 @@
 # Progression of infection within individuals
 #
 import random
+import typing
+
 import numpy as np
 from collections import defaultdict
 
@@ -25,9 +27,11 @@ def __init__(self):
         to a current infection status of a person. The columns of that
         row then indicate the transition probabilities to the remaining
         infection statuses. Number of infection states is set by
-        taking the size of the InfectionStatus enum. Transition time matrix
-        is also initialised and associated parameters are called from the
-        parameters class.
+        taking the size of the InfectionStatus enum. Waning transition matrix
+        is initialised and adapts the state transition matrix with rate
+        multipliers relating to waning immunity which are called from the
+        parameters class. Transition time matrix is also initialised and
+        associated parameters are called from the parameters class.
 
         Infectiousness progression defines an array used to scale a person's
         infectiousness and which depends on time since the start of the
@@ -38,8 +42,13 @@ def __init__(self):
         use_ages = Parameters.instance().use_ages
         coefficients = defaultdict(int, Parameters.instance()
                                    .host_progression_lists)
-        matrix_object = StateTransitionMatrix(coefficients, use_ages)
+        multipliers = defaultdict(list, Parameters.instance()
+                                  .rate_multiplier_params)
+        matrix_object = StateTransitionMatrix(coefficients, multipliers,
+                                              use_ages)
         self.state_transition_matrix = matrix_object.matrix
+        if pe.Parameters.instance().use_waning_immunity:
+            self.waning_transition_matrix = matrix_object.waning_matrix
 
         self.number_of_states = len(InfectionStatus)
         assert self.state_transition_matrix.shape == \
@@ -132,7 +141,7 @@ def set_infectiousness(person: Person, time: float):
         if person.infection_start_time < 0:
             raise ValueError('The infection start time cannot be negative')
 
-    def update_next_infection_status(self, person: Person):
+    def update_next_infection_status(self, person: Person, time: float = None):
         """Assigns next infection status based on current infection status
         and on probabilities of transition to different statuses. Weights
         are taken from row in state transition matrix that corresponds to
@@ -143,8 +152,11 @@ def update_next_infection_status(self, person: Person):
 
         Parameters
         ----------
-        Person : Person
+        person : Person
             Instance of person class with infection status attributes
+        time : float
+            Current simulation time (if necessary for the method, default =
+            None)
 
         """
         if person.infection_status in [InfectionStatus.Dead,
@@ -167,11 +179,24 @@ def update_next_infection_status(self, person: Person):
                 person.next_infection_status = InfectionStatus.Dead
                 return
         row_index = person.infection_status.name
-        weights = self.state_transition_matrix.loc[row_index].to_numpy()
-        weights = [w[person.age_group] if isinstance(w, list) else w
-                   for w in weights]
-        outcomes = range(1, self.number_of_states + 1)
 
+        # If we are not using waning immunity or person.time_of_recovery is
+        # None (so they have never reached Recovered) then we choose weights
+        # from the state_transition_matrix. Otherwise, we use the
+        # waning_transition_matrix.
+        if (not Parameters.instance().use_waning_immunity or
+                not person.time_of_recovery):
+            weights = self.state_transition_matrix.loc[row_index].to_numpy()
+            weights = [w[person.age_group] if isinstance(w, list) else w
+                       for w in weights]
+        else:
+            if time is None:
+                raise ValueError("Simulation time must be passed to "
+                                 "update_next_infection_status when waning "
+                                 "immunity is active")
+            weights = self._get_waning_weights(person, time)
+
+        outcomes = range(1, self.number_of_states + 1)
         if len(weights) != len(outcomes):
             raise AssertionError('The number of infection statuses must' +
                                  ' match the number of transition' +
@@ -183,6 +208,46 @@ def update_next_infection_status(self, person: Person):
 
         person.next_infection_status = next_infection_status
 
+    def _get_waning_weights(self, person: Person, time: float) -> typing.List:
+        """Given that the current person has previously recovered, this method
+        will return a list of updated weights based on the level of immunity
+        the person has. The weights taken from the waning_transition_matrix
+        are lambda expressions parameterized by t (time_since_recovery).
+
+        Parameters
+        ----------
+        person : Person
+            Instance of person class with infection status attributes
+        time : float
+            Current simulation time
+
+        Returns
+        -------
+        list:
+            List of weights representing the probability of transitioning to
+            a given compartment.
+        """
+        row_index = person.infection_status.name
+        time_since_recovery = time - person.time_of_recovery
+        weights = list(self.waning_transition_matrix.loc[row_index])
+
+        # Note that below, each entry w will be a lambda expression returning
+        # a np.array either representing a float (shape = ()) or a list
+        # (shape = (n,)) hence the conditions.
+        new_weights = []
+        for w in weights:
+            if isinstance(w, (int, float)):
+                new_weights.append(w)
+            else:
+                # This is evaluating the lambda expressions at t =
+                # time_since_recovery
+                w_evaluated = w(time_since_recovery)
+                if w_evaluated.shape:
+                    new_weights.append(w_evaluated[person.age_group])
+                else:
+                    new_weights.append(w_evaluated)
+        return new_weights
+
     def update_time_status_change(self, person: Person, time: float):
         """Calculates transition time as calculated in CovidSim,
         and updates the time_of_status_change for the given
@@ -193,7 +258,7 @@ def update_time_status_change(self, person: Person, time: float):
 
         Parameters
         ----------
-        Person : Person
+        person : Person
             Instance of Person class with :class:`InfectionStatus` attributes
         time : float
             Current simulation time
@@ -278,13 +343,12 @@ def _updates_infectiousness(self, person: Person, time: float):
             person.infectiousness = person.initial_infectiousness *\
                 scale_infectiousness[time_since_infection]
         # Sets infectiousness to 0 if person just became Recovered, Dead, or
-        # Vaccinated, and sets its infection start time to None again.
+        # Vaccinated
         elif person.infectiousness != 0:
             if person.infection_status in [InfectionStatus.Recovered,
                                            InfectionStatus.Dead,
                                            InfectionStatus.Vaccinated]:
                 person.infectiousness = 0
-                person.infection_start_time = None
 
     def __call__(self, time: float):
         """Sweeps through all people in the population, updates their
@@ -319,7 +383,7 @@ def __call__(self, time: float):
                         self.set_infectiousness(person, time)
                         if not person.is_symptomatic():
                             asympt_or_uninf_people.append((cell, person))
-                    self.update_next_infection_status(person)
+                    self.update_next_infection_status(person, time)
                     if person.infection_status == InfectionStatus.Susceptible:
                         person.time_of_status_change = None
                         break

pyEpiabm/pyEpiabm/sweep/initial_infected_sweep.py

@@ -104,6 +104,7 @@ def __call__(self, sim_params: dict):
                                            ["initial_infected_number"]))
         for person in pers_to_infect:
             person.update_status(InfectionStatus.InfectMild)
+            person.increment_num_times_infected()
             person.household.remove_susceptible_person(person)
             person.next_infection_status = InfectionStatus.Recovered
             HostProgressionSweep.set_infectiousness(person, start_time)