In order to capture the interactive dynamics that emerge with more than one cell type available, the single cell population differential equation model with delayed viral production rates, as it has been proposed in (Baccam 2006), was extented to a two target cell model in (Dobrovolny 2010). In our models, we additionally allowed each virus entity to mutate mutate from its wildtype into a more drug-resistant instance. In terms of the nomenclature, we separated the two target cells into default (subscript d) and secondary (subscript s) cells with each containing a wildtype (subscript wt) and a mutant (subscript \(\mu\)) sub-population.

Initial values of wildtype virus (\(V_{\text{wt}}\)) and mutant virus (\(V_{\mu}\)) proceed to infect primary target cells (\(T_{\text{d}}\)) at rate \(\beta\) and secondary target cells (\(T_{\text{s}}\)) at rate \(r_{\beta}\beta\). Once infected, cells migrate into their eclipse phase (\(E\)) and then turn into productively infected cells \(I\) at rates \(\tau_{\text{E}}^{-1}\) and \(\tau_{\text{I}}^{-1}\), respectively. Of both the eclipse and the infected cells there are four distinct types: any combination of default and secondary with wildtype and mutant represent possible states. Once primary (secondary) target cells have reached their productive stage, they will produce virus at rate \(p\) (\(r_{p}p\)) while slowly dying off at rate \(c\). When dying they accumulate as \(D\), from which they may regenerate back to available target cells \(T\) at rate \(\ell\).
The resulting model is a system comprised of 14 differential equations, see eq. (\ref{eq:ode_system_T}) to (\ref{eq:ode_system_D}).

\begin{align}
\label{eq:ode_system_D}
\dot{\mathbf{T}}= & \,\begin{pmatrix}{}\dot{T}_{\text{d}}\\
\dot{T}_{\text{s}}\end{pmatrix}=-\big{[}\beta_{\text{wt}}V_{\text{wt}}+\beta_{\mu}V_{\mu}\big{]}\begin{pmatrix}1&0\\
0&r_{\beta}\end{pmatrix}\mathbf{T}+\ell\mathbf{D} \\
\dot{\mathbf{E}}_{j}{}= & \,\begin{pmatrix}\dot{E}_{\text{d}}^{j}\\
\dot{E}_{\text{s}}^{j}\end{pmatrix}=(1-m_{j})\beta_{j}V_{j}\begin{pmatrix}1&0\\
0&r_{\beta}\end{pmatrix}\mathbf{T}-\tau_{E}^{-1}\mathbf{E}_{j}\hskip 20.485984pt\text{where}\hskip 6.259606ptj=\text{wt}\,,\mu \\
\dot{\mathbf{I}}_{j}{}= & \,\begin{pmatrix}\dot{I}_{\text{d}}^{j}\\
\dot{I}_{\text{s}}^{j}\end{pmatrix}=\tau_{E}^{-1}\mathbf{E}_{j}-\tau_{I}^{-1}\mathbf{I}_{j} \\
\dot{V}_{j}{}= & \,(1-n_{j})p_{j}\begin{pmatrix}1-\mu_{\text{nt}}\\
r_{p}-r_{p}\mu_{\text{nt}}\\
\mu_{\text{nt}}\\
r_{p}\mu_{\text{nt}}\end{pmatrix}^{\intercal}\begin{pmatrix}\mathbf{I}_{\text{wt}}\\
\mathbf{I}_{\mu}\end{pmatrix}-cV_{j} \\
\label{eq:ode_system_D}\dot{\mathbf{D}}{}= & \,\begin{pmatrix}\dot{D}_{\text{d}}\\
\dot{D}_{\text{s}}\end{pmatrix}=\tau^{-1}\big{(}\mathbf{I}_{\text{wt}}+\mathbf{I}_{\mu}\big{)}-\ell\mathbf{D}\\
\end{align}

The key parameters that control the differences between the two cell populations are the relative susceptibility to infection (\(r_{\beta}\in\mathbb{R}^{+}\)), the relative viral production rate (\(r_{p}\in\mathbb{R}^{+}\)) and the fraction of initial secondary target cells (\(r_{T}\in[0,1]\)). Note that \(r_{T}\) is only used in setting up the initial target cell population and does not explicitly appear in the system of differential equations.

Throughout simulations, the initial number of overall target cells was set to \(T_{0}=4e8\), which coincides with anatomical estimates for the human upper respiratory tract (Baccam 2006). Accordingly, we fixed \(r_{\text{T}}\) at \(70\%\) equating to the proportion of default cells that predominantly express \(\text{SA}\alpha\text{2,6}\,\text{Gal}\) receptors on their surface. This is based on a study of human lung physiology indicating that the epithelium of the upper airway (up to the fifth generation) comprises 50–85% nonciliated cells (Crystal 1991), and reports indicating that nonciliated cells predominantly express \(\text{SA}\alpha\text{2,6}\,\text{Gal}\) receptors on their surface (Thompson 2006, Ibricevic 2006, Kogure 2006). The remaining \(30\%\) of secondary cells are identified to express \(\text{SA}\alpha\text{2,6}\,\text{Gal}\) receptors on their surface. More of the default parameter choices that were made in the simulations are listed in Table \ref{table:parameters}.

Symbol | Description | Default value |
---|---|---|

\(T_{0}\) | Total number of target cells | \(4e8\) |

\((E_{0},I_{0},D_{0})\) | Initially eclipsed, infected, dead cells | \((-,-,-)\) |

\(V_{0}\) | Initial viral inoculum | \(7.5e-2[V]\) |

\(\tau_{E}\) | Length of eclipse phase | \(6h\) |

\(\tau_{I}\) | Length of virus production phase | \(\approx 4.36h\) |

\(c\) | Virus clearance rate | \(\approx 4.36h\) |

\(\beta_{\text{d,wt}}\) | Cell infection rate (default, wildtype) | \(3.2e-5[V]^{-1}\cdot d^{-1}\) |

\(p_{\text{d,wt}}\) | Virus production rate (default, wildtype) | \(4.6e-2[V]^{-1}\cdot d^{-1}\) |

Our two cell cell population model is able to capture a large variety of cell receptor specifics and, in other words, is thus situated in a large parameter space. In a recent work, this space has been studied in detail for the specific parameters that are able to induce a severe influenza infection (Dobrovolny 2010). Since it is the purpose of this study to assess the efficacy of drugs when applied to such a long lasting infection, and in order to enforce a severe infection (lasting for significantly more than five days), we chose (unless otherwise noted) the ratio between the default and secondary cell population’s production rate (\(r_{\beta}\)) and the value that indicating the relative virus release rate of the default and the secondary population (\(r_{p}\)) as \((r_{\beta},r_{p})=(1e-4,2e3)\). For a comparison of the resulting severe influenza with a common influenza infection, see Fig. \ref{fig:std_vs_severe}.

In order to study the likelihood of an infection impairing its host as a function of the fitness of the mutant virus population, we measured the probability of the viral titer amount surpassing a population value of 1e8, see Fig. \ref{fig:breakthrough_infections}. As illustrated in Fig. \ref{fig:ode_system}, an adamantanes (ADM) based drug prevents the uncoating of the virus, and in order to study its efficacy in different infections we varied the relative ratio \(\beta_{\text{fit}}\) between mutant and wildtype infection rates. When studying the efficacy of a neuraminidase inhibitors (NAIs) based drug (that blocks the release of viral particles), we varied instead how the release rate of mutant and wildtype virus compare (denoted by \(p_{\text{fit}}\)). Both measurements were repeated which were able to identify a similar general trend, regardless of the precise potency of the drug.

\begin{equation}
\beta_{\text{fit}}=\frac{\beta_{\mu}}{\beta_{\text{wt}}}\hskip 28.452756ptp_{\text{fit}}=\frac{p_{\mu}}{p_{\text{wt}}}\\
\end{equation}

\label{fig:breakthrough_infections}{subfigure}
[b]0.45 {subfigure}[b]0.45

Figure 1: The first plot shows how likely an adamantane based drug is to supress an infection. We varied with what rate the mutant virus will infect available target cells, compared to the wildtype virus. Also, the efficacy \(\varepsilon_{\text{ADM}}\) was varied to confirm the overall trend of a higher likelihood of breakthrough infections when the mutant virus populations infection rate is increased. The second plot similarly shows the supressing behavior of an NAI based drug for varied virus release rates in mutant and wildtype virus populations. Again, simulations were repeated for different drug efficacies \(\varepsilon_{\text{NAI}}\).To evaluate the average contribution of the mutant virus to an overall infection, we monitored the fraction of mutant virus produced, relative to the total virus amount, see Fig. \ref{fig:mutant_virus_amount}. Simulations were repeated for two choices of parameters, depending on the drug type that was being applied:

- •
For an adamantane based drug, we varied the infection rate by tuning \(p_{\text{fit}}\).

- •
We varied the viral production rate when studying NAI based drugs.

We expect the virus levels in an infection to surpass detection levels (which we fixed at a total number of \(10^{5}\) viral particles) quicker, if the mutant virus were to be fitter than the wildtype virus. We studied this by monitoring the times of detection (\(t_{\text{det}}\)) in dependency of either \(\beta_{\text{fit}}\) (in case of adamantane based drugs) or \(p_{\text{fit}}\) (when administering times of detection in an infection challenged by NAI based drugs). From Fig. \ref{fig:times_of_detection} we can observe that higher efficacies result in a suppression of the development of the infection, making times of detection become greater. In addition, when the efficacy approaches zero, the measured times of detection become centered around its mean value.

[b]0.45 {subfigure}[b]0.45

Figure 3: In the first figure, times of detection \(t_{\text{det}}\) are plotted against changing values of \(\beta_{\text{fit}}\) while an adamantane based treatment is in effect. In the second figure, similarly, the times of detection are evaluated for a given \(p_{\text{fit}}\), while an NAI based drug is being applied. The shaded area that surround both plots indicate the standard deviation of each of the individual measurements and the colors signal the corresponding efficacy. The black curve in both plots shows that times of detection stabilize when setting the efficacy to zero.To optimize the dampening effect of adamantane or NAI based drugs on the evolution of virus levels in an influenza infection we tested a set of periodic treatments. We control the treatment using two parameters: The drug is first being applied for a time \(t_{\text{on}}\). Next, the treatment is paused for a time \(t_{\text{off}}\). The treatment is then periodically turned on and off again. As an example of how infection curves behave for different parameter choices, see Fig. \ref{fig:periodic_treatment}.

[b]0.3
{subfigure}[b]0.3
{subfigure}[b]0.3

{subfigure}[b]0.3
{subfigure}[b]0.3
{subfigure}[b]0.3

We now varied the parameter that controls the regeneration of dead cells, \(\ell\). Once they come back to life, they reappear as available target cells. For different choices of the regeneration rates, we now monitored the total infection time, ie. the total time the virus population was larger than the breakthrough infection level of 10e8 viral particles (cf. Fig. \ref{fig:breakthrough_infections}). All measurements were recorded in the severe infection regime (Dobrovolny 2010), thus our parameter choices for the fitness ratios between default and secondary populations were fixed at \((r_{\beta},r_{p})=(10^{-4},10^{3.5})\). The relative amount of initial mutant virus was set to \(r_{\text{T}}=0.3\) throughout simulations.

The different periodic treatment schemes, which get represented by the pairings \((t_{\text{on}},t_{\text{off}})\), span the parameter space we are interested in. In order to account for a variety of schemes, we chose all pairings between \(1dpi\) and \(10dpi\), with a stepsize of \(1/4dpi\). Again, for some exemplary measurements see Fig.

[b]0.28
{subfigure}[b]0.28
{subfigure}[b]0.28

{subfigure}[b]0.28
{subfigure}[b]0.28
{subfigure}[b]0.28

In order to measure how much sense it makes to vary the treatment scheme in periodic drug applications, we looked at how strongly the total infection time varies for fixed \(t_{\text{on}}\), when varying \(t_{\text{off}}\). We then averaged over each of these observations to arrive at the observed coefficient of variation as a function of the regeneration rate, \(\ell\).

- •
We first calculate the average total infection time for fixed \(t_{\text{on}}\)’s and denote that by \(\overline{\text{TIT}}_{\ell}(t_{\text{on}})\).

- •
Also, we are interested in the variance of the above value, when varying \(t_{\text{off}}\).

\begin{equation} \operatorname{Var}[\text{TIT}_{\ell}(t_{\text{on}})]=\frac{1}{N}\sum_{\{t_{\text{off}}\}}\Big{(}\text{TIT}_{\ell}(t_{\text{on}},t_{\text{off}})-\overline{\text{TIT}}_{\ell}(t_{\text{on}})\Big{)}^{2}\nonumber \\ \end{equation} - •
The two values are combined into the coefficient of variation for a given \(\ell\) and \(t_{\text{on}}\).

\begin{equation} c_{\text{v},\ell}(t_{\text{on}})=\frac{\sqrt{\operatorname{Var}[\text{TIT}_{\ell}(t_{\text{on}})]}}{\overline{\text{TIT}}_{\ell}(t_{\text{on}})}\nonumber \\ \end{equation} - •
The coefficient of variation is a natural measure of the amount of influence that the variation of \(t_{\text{off}}\) has. Thus, the larger \(c_{\text{v},\ell}(t_{\text{on}})\), the more susceptible we assume the infection to be to variations in the periodic treatment scheme.

- •
In order to study the coefficient of variation as a function of \(\ell\) only, one can average over all \(t_{\text{on}}\). This yields:

\begin{equation} \overline{c}_{\text{v},\ell}=\frac{1}{N}\sum_{\{t_{\text{on}}\}}c_{\text{v},\ell}(t_{\text{on}})\nonumber \\ \end{equation} - •
In the end, a large coefficient of variation can be identified with a large influence of periodic treatment schemes. If \(\overline{c}_{\text{v},\ell}\) is small, periodic treatments have little effect.