Measles-ish SEIRS Models in Python

Model 1

We're going to explore how for an ODE-type model, with a constant birth rate (supply of new susceptibles), and initial conditions that are endemic measles-ish, we get a damped oscillator as output.

The periodicity and rate of damping would both be interesting to explore.

Cf. http://math.uchicago.edu/~shmuel/Modeling/Keeling%20and%20Rohani/chap%202.pdf. "An important issue for any dynamical system concerns the manner in which a stable equilibrium is eventually approached. Do trajectories undergo oscillations as they approach the equilibrium state or do they tend to reach the steady state smoothly? The SIR system is an excellent example of a 'damped oscillator,' which means the inherent dynamics contain a strong oscillatory component, but the amplitude of these fluctuations declines over time as the system equilibrates (Figure 2.5). Figure 2.6, shows how the period of oscillations (as determined by equation (2.23) in Box 2.4) changes with the transmission rate (β) and the infectious period (1/γ ). We note that (relative to the infectious period) the period of oscillations becomes longer as the reproductive ratio approaches one; this is also associated with a slower convergence toward the equilibrium." (Figures not shown)

Crude Birth Rate

42 / 50

beta

2 / 10

import matplotlib.pyplot as plt # Parameters population_size = 1e6 initial_infected = 1000 #initial_susceptible = 999000 beta = BETA # Transmission rate gamma_e = 1/7 # Rate of transition from exposed to infectious gamma_i = 1/7 # Rate of transition from infectious to recovered cbr = CBR # Calculate births per day births_per_day = ((population_size / 1000) * (cbr / 365)) #KM: Just a note that this works with constant population, but if we want to accurately capture compounding, #then this will accrue small errors over time. Not an issue here, just a note that in general, daily CBR = (1+CBR)^(1/365) - 1 # Initial state E = 0 I = initial_infected R = 0.9*population_size S = population_size-R-I #KM: As you sweep over beta, you could start closer to equilibrium and speed up burnin, #by using R = (1-1/R0)*pop, or ~ (1-gamma/beta)*pop. # Simulation parameters total_time = 365*20 # Total simulation time in days timestep = 1 # Time step size in days # Lists to store compartment counts over time susceptible_over_time = [S] exposed_over_time = [E] infectious_over_time = [I] recovered_over_time = [R] time_steps = [0] # Simulation loop for t in range(total_time): # Process births births = births_per_day # Calculate transitions new_exposed = beta * S * I / population_size new_infectious = E * gamma_e new_recovered = I * gamma_i # Update state S += births - new_exposed E += new_exposed - new_infectious I += new_infectious - new_recovered R += new_recovered #KM: Total population is growing, but I don't think it matters. # Ensure non-negative population counts S = max(S, 0) E = max(E, 0) I = max(I, 0) R = max(R, 0) # Store compartment counts at current timestep # KM: My experience in Matlab would lead me to believe that pre-allocating these vectors might give a substantial speedup susceptible_over_time.append(S) exposed_over_time.append(E) infectious_over_time.append(I) recovered_over_time.append(R) time_steps.append(t+1) # Check for end of epidemic if I == 0: print("Epidemic ended.") break # Plot compartment counts over time fig, ax1 = plt.subplots(figsize=(10, 6)) # Plot susceptible, exposed, and infectious on the primary y-axis ax1.plot(time_steps, susceptible_over_time, label='Susceptible', color='blue') ax1.plot(time_steps, exposed_over_time, label='Exposed', color='orange') ax1.plot(time_steps, infectious_over_time, label='Infectious', color='red') ax1.set_xlabel('Time (days)') ax1.set_ylabel('Population count') ax1.set_title('SEIR Model Simulation') # Create a second y-axis for recovered ax2 = ax1.twinx() ax2.plot(time_steps, recovered_over_time, label='Recovered', color='green') ax2.set_ylabel('Recovered count') # Show legend for both axes ax1.legend(loc='upper left') ax2.legend(loc='upper right') plt.title('SEIR Model Simulation') plt.legend() plt.grid(True) plt.show()

Run to view results

# Let's draw S-I orbitals aka a scattler plot of S vs I import numpy as np plt.scatter(susceptible_over_time, np.add( infectious_over_time, exposed_over_time )) plt.title("S-I orbital") plt.xlabel("S") plt.ylabel("I + E") plt.show()

Run to view results

You can easily play with different values of beta & birth rate and see how it varies the periodicity of the infectious curve and the rate at which it "damps". But the main takeaway here is that this type of model converges at a steady-state (in the time domain).

Questions

What determines the periodicity? (The formula for this is 2*pi*sqrt(average_age_at_infection * generation_time)

What determines the amplitude? Damping rate?

Does this reproduce in ABMs? (See below)

What is the role of the fact that this model is all continuous values for the number of people? (It means infection can still survive at tiny levels because really those values would round to zero infected individuals in practice.)

What is the role of the lack of stochasticity/noise? (Cf. Keeling & Rohani Ch. 6: "As most SIR-type disease models approach equilibrium in a series of decaying epidemics, the increased transient-like dynamics of stochastic models leads to oscillations close to the natural frequency (Chapter 2). Thus, stochasticity can excite epidemic oscillations around the normal endemic prevalence, leading to sustained cycles.")

Run to view results

Model 2. CCS Quick Explorer

Here we create a toy SEIRS Endemic Measles-ish model with discrete number of individuals. We stop early upon elimination. Initially this was deterministic but we've added some stochasticity to avoid unhelpful artifacts that can emerge from a purely deterministic approach.

beta

1 / 20

cbr

10 / 50

Initial Population

1000000 / 1000000

import numpy as np import matplotlib.pyplot as plt import sys from collections import deque # Parameters population = POP # Total population cbr = CBR # Crude birth rate per 1000 people per year beta = BETA # Transmission rate incubation_period = 7 # Incubation period in days infectious_period = 7 # Infectious period in days simulation_days = 7300*3 # Number of days to simulate # Calculate births per day births_per_day = int((population / 1000) * (cbr / 365)) # Initial values E_queue = deque([1000] * incubation_period) # Exposed individuals queue I_queue = deque([0] * infectious_period) # Infectious individuals queue (seeded outbreak) S = [int(population/(beta*infectious_period))] R = [population - S[0] - sum(E_queue) - sum(I_queue)] # Make reporting lists for E and I since we're modeling with queues E_reporting = [sum(E_queue )] I_reporting = [sum(I_queue )] # Simulation for day in range(1, simulation_days): # Calculate new exposures # We can cancel a couple of populations from this math but that's not necessarily the main point yet. old_exposures = int((S[-1] / population) * (sum(I_queue) / population) * beta * population) #KM: where does old exposures get used? new_exposures = np.round(beta * sum(I_queue) * S[-1] / population) new_exposures = np.random.poisson(new_exposures) # Calculate new infections (from the exposed population) new_infections = E_queue.pop() # Push new exposures into the exposed queue E_queue.appendleft(new_exposures) # Push new infections into the infectious queue I_queue.appendleft(new_infections) # Update recovered population R.append(R[-1] + I_queue.pop()) #KM: I'm not 100% familiar with queues, but does it matter that we are popping off the E_queue before #appending new exposures, but appending new infections before popping it? I don't think so but curious. # Update susceptible population S.append(S[-1] - new_exposures + births_per_day) population += births_per_day E_reporting.append( sum(E_queue ) ) I_reporting.append( sum(I_queue ) ) if sum(I_queue)+sum(E_queue) == 0: print( "Eliminated!" ) break time_steps = np.arange(len(S)) fig, ax1 = plt.subplots(figsize=(10, 6)) # Plot susceptible, exposed, and infectious on the primary y-axis ax1.plot(time_steps, E_reporting, label='Exposed', color='orange') ax1.plot(time_steps, I_reporting, label='Infectious', color='red') ax1.set_xlabel('Time (days)') ax1.set_ylabel('Population count') # Create a second y-axis for recovered ax2 = ax1.twinx() ax2.plot(time_steps, S, label='Susceptible', color='blue') ax2.plot(time_steps, R, label='Recovered', color='green') ax2.set_ylabel('Recovered count') # Combine legends from both axes lines1, labels1 = ax1.get_legend_handles_labels() lines2, labels2 = ax2.get_legend_handles_labels() lines = lines1 + lines2 labels = labels1 + labels2 ax1.legend(lines, labels, loc='upper right') plt.xlabel('Days') plt.ylabel('Population') plt.title('SEIR Model Simulation') #KM: Might I suggest: plt.title('SEIR Model Simulation: beta = '+str(beta)+', CBR = '+str(cbr)+', pop = '+str(population)) #That way the plot shows the simulation parameters plt.show()

Run to view results

We expect to see that for sub-CCS population sizes we eliminate. For others we reach endemecity.

Depending on some of the values in our model, the epidemic will 'burn itself out' (eliminate) or persist. The most important values which determine that are beta, cbr, and population. For a given beta & cbr, there will be a threshold population. For example, with beta=2 and cbr=25, the threshold population for CCS appears experimentally to be about 375,000.

You can adjust those values yourself and see what you find.

Note that there is no stochasticity in the above model, which means we can repeatedly maintain endemecity at lower populations by going really "close to the line" during the low prevalence periods. Adding stochasticity will result in more frequent elimination, meaning the true CCS population threshold will be higher than what this toy example shows.

CCS Sweeps & Surface

Now let's run a sweep across a range of beta and cbr values of interest and see at what population we hit CCS.

Run to view results

Note that these values seem higher (for low CBR settings) than I've heard before.

.css-15w88e5{color:var(--chakra-colors-fg-neutral-primary);font-weight:inherit;letter-spacing:-0.09px;}Measles-ish SEIRS Models in Python

Model 1

Questions

Model 2. CCS Quick Explorer

CCS Sweeps & Surface

TBD: Impact of Seasonality on CCS

Measles-ish SEIRS Models in Python