# Exercise: Modeling multiple strains

Suppose we have time series of cases caused by two strains, and we want to know if the strains interact. To develop expectations of the possible behaviors of two strains, we build a simple, ODE-based model. The one below has been adapted from Keeling & Rohani (2007). Subscript indices refer to the status with respect to strains 1 and 2, respectively.

(1) Draw the compartmental diagram for this model, adding parameters to the flows. Interpret each parameter. What kinds of restrictions (e.g., bounded between 0 and 1) make sense for each? Is this a history-based or status-based model?

(2) What is the model assuming about the duration and strength of interaction? What other potential modes of interaction are missing from this model?

(3) Consider two strains of your favorite antigenically variable pathogen. If you wanted to make a model that captures the basic dynamics of these strains, which assumptions of this model would you relax first?

(Optional) Solve analytically for the equilibria and conditions for invasion for each strain when the other is at its endemic equilibrium.

Compare answers with your neighbors, and be sure to raise any interesting points of disagreement (or agreement or confusion) with the class.

For this next section you’ll need to run a Python script.

This can either be downloaded here and run locally or can be run directly from within MyBinder at . This script outputs the figure `fig_plot_two_strain_ts.png`

.

Let’s assume we think this model is a pretty good starting point for testing hypotheses.
We’ll numerically integrate these equations to see what they predict.
Open `two_strain.py`

and `plot_two_strain_ts.py`

and examine the code.
It uses the Euler method to solve the system.
The function calculates the state variables over time at `step_size`

intervals.
Note that it produces estimates of the state variables over a different schedule (`t=0`

to `t=end_time`

at intervals of `output_interval`

).
In addition to the step sizes for integration and outputs, the `two_strain`

function requires parameters and initial conditions.

Run this with:

```
python plot_two_strain_ts.py
```

(4) Modify the parameters so that the strains do not interact, and examine the resulting time series. How can you confirm that the integration is accurate?

(5) Without allowing for interaction, see if you can adjust the parameters to obtain the following scenarios: (i) the strains have biennial, out-of-phase dynamics, and (ii) both strains are chaotic. Do you see any recurring patterns in the chaotic dynamics?

(6) Now change

`alpha`

and`a`

to make the strains interact. How do the dynamics differ? Do you see differences in the dynamics when you change`alpha`

v.`a`

?

(7) Let’s use the model to consider the conditions under which a new variant can invade—specifically, the invasion of Omicron in late 2021. (Pretend there was a single resident strain, Delta, and a single invader, Omicron.) Major uncertainties at the time included the fraction of the population immune to infection with Omicron compared to prior variants such as Delta, the intrinsic reproductive numbers of Omicron and Delta, and the cross-immunity after infection with each to the other. In laboratory experiments, sera from recipients of the mRNA-1273 vaccine (containing the original Wuhan strain, estimated to be half as transmissible as Delta) experienced a 30-fold reduction in neutralization against Omicron (the reduction was ~2.5-fold for Delta). Epidemiological studies in England suggested that Omicron was 50\% more transmissible than Delta within households and potentially twice as transmissible outside households. By adjusting the parameters in

`two_strain.py`

, what intrinsic reproductive numbers, differences in immunity to Delta vs. Omicron (i.e., degree of immune escape in Omicron), and values of`a`

or`alpha`

are consistent with Omicron’s rapid displacement of Delta? Try to identify a few contrasting scenarios.

(Optional) Consult the documentation for SciPy’s’ ODE solver

`odeint`

in`scipy.integrate`

. This solver uses an adaptive step size. See if you can implement the model using the solver–most of the content of`two_strain.py`

can be copied and pasted. Investigate the consequences of changing the tolerance limits (`rtol`

and`atol`

) and acceptable step sizes (`hmin`

and`hmax`

) of the solver. These equations are said to be stiff.

As you’ve probably discovered, SIR-type models with seasonal forcing can generate diverse patterns.
To summarize these patterns, we’ll create bifurcation diagrams over several of the parameters we’re interested in.
Bifurcation diagrams show the *equilibrium* behavior of a system as one parameter is varied.
These summaries are often “strobes” of the system (e.g., the number of cases every Jan. 1), but they can also show all maxima and minima over some period.
Often, for each value of the varied parameter, the system is simulated to equilibrium from different starting conditions.
This approach can reveal the presence of multiple stable states or attractors, as shown in the diagram below for a single-strain model (Earn et al. 2001).

For simplicity, we’ll create our bifurcation diagrams using the same initial conditions throughout. We won’t measure precisely whether we have obtained equilibrium conditions–we’ll simply simulate for a long time, and sample the few final time points.

(8) Carefully choose reasonable default values for all your parameters. Then create bifurcation diagrams with

`plot_two_strain_bd.py`

, independently sweeping over`alpha`

,`a`

, and`beta`

for one strain. Be modest in creating your first diagram, and record how long it takes to run. Check your results by simulating a few time series at different points.

Run this with:

```
python plot_two_strain_bd.py
```

(9) Is reduced transmission following infection with the heterologous strain dynamically equivalent to reduced susceptibility?

(10) If you’re starting with seasonal forcing, what happens when you remove it? Do the reverse exercise if you start without it.

You might be wondering why we should bother with all this exploration. (Can’t we just fit the darn things already?) Fitting fails all the time, and it’s important to understand why. Simulating and summarizing dynamics are also useful tools to investigate the results from inference.

(11) Let’s say you obtained an estimate of

`beta`

from fitting the model. You make a bifurcation diagram, and it indicates that a 5% increase in`beta`

would shift the system to chaotic dynamics. The time series don’t show chaotic dynamics. What do you conclude? What if your estimate of`beta`

came from experiment?