Modeling the measles paradox reveals the importance of cellular immunity in regulating viral clearance

Measles virus (MV) is a highly contagious member of the Morbillivirus genus that remains a major cause of childhood mortality worldwide. Although infection induces a strong MV-specific immune response that clears viral load and confers lifelong immunity, transient immunosuppression can also occur, leaving the host vulnerable to colonization from secondary pathogens. This apparent contradiction of viral clearance in the face of immunosuppression underlies what is often referred to as the ‘measles paradox’, and remains poorly understood. To explore the mechanistic basis underlying the measles paradox, and identify key factors driving viral clearance, we return to a previously published dataset of MV infection in rhesus macaques. These data include virological and immunological information that enable us to fit a mathematical model describing how the virus interacts with the host immune system. In particular, our model incorporates target cell depletion through infection of host immune cells—a hallmark of MV pathology that has been neglected from previous models. We find the model captures the data well, and that both target cell depletion and immune activation are required to explain the overall dynamics. Furthermore, by simulating conditions of increased target cell availability and suppressed cellular immunity, we show that the latter causes greater increases in viral load and delays to MV clearance. Overall, this signals a more dominant role for cellular immunity in resolving acute MV infection. Interestingly, we find contrasting dynamics dominated by target cell depletion when viral fitness is increased. This may have wider implications for animal morbilliviruses, such as canine distemper virus (CDV), that cause fatal target cell depletion in their natural hosts. To our knowledge this work represents the first fully calibrated within-host model of MV dynamics and, more broadly, provides a new platform from which to explore the complex mechanisms underlying Morbillivirus infection.


Viral load
As measles virus (MV) infection causes a cell-associated viremia, viral load was determined as the number of productively infected cells per million peripheral blood mononuclear cells (PBMC). Using an infectious center assay, several replicates of multiple tenfold dilutions of PBMC were tested for their ability to infect susceptible cells. Cocultivation was screened between 4 and 7 days and a well was either scored as positive or negative for infection. The 50% point was then determined using the formula of Reed and Muench [2], and the result was normalized to the number of infected cells per million PBMC based on the number of cells used in the first dilutions.

Lymphocytes
We use the total number of lymphocytes as a proxy for the main virus target cells. Lymphocyte counts were always collected in the morning, and the total number of lymphocytes per microliter of blood was measured using an automated cell counter. NK and NKT cells are not (or at best poorly) susceptible to MV infection and are not expected to contribute substantially to this population. In addition, monocytes and granulocytes were distinguished as separate populations. Thus the majority of the total lymphocyte counts in our data represent cells susceptible to infection.

MV-specific T cells
The IFN-γ spot forming cells (SFCs) were determined after MV-specific stimulation, and data shown are specific counts after the subtraction of counts obtained following mock stimulation. The data therefore represent antigen-specific cells and are used as a proxy for MV-specific activated T cells.

Avidity
Avidity was measured by enzyme immunoassay using plates coated with lysates from MeV-infected cells. An avidity index was calculated by the ability of MeV-specific plasma IgG to be eluted by increasing concentrations of NH4SCN (see for example [3]).

Alternative model formulations T cell activation
The timing and magnitude of T cell activation is modeled in the main text using a continuous saturation function described previously for LCMV [4]. However, we also investigated alternative functional forms. First, we defined activation within a discrete timeframe such that with all other model equations remaining unchanged. Note that this function has two parameters that must be estimated (t 1 and t 2 ), compared to one parameter for the continuous saturation function (s). Second, we investigated activation that was proportional to viral load, i.e. we set the function governing T cell proliferation in Eq 3 to f (V ) = V , and all other terms containing 1 − f (V ) to 1. This formulation has no additional parameters to estimate and results in the following system of equations:

Lymphocyte proliferation
Proliferation of general susceptible lymphocytes (S) is modeled in the main text using a discrete on/off function to describe early (but temporary) lymphocyte activation i.e. q sδ (t)S, whereδ Here, proliferation is controlled by two parameters: the rate of proliferation, q s , and the duration of proliferation, t d . However, we also investigated two alternative functional forms with fewer parameters. First, we defined proliferation to be constant and solely dependent on the parameter q s by settingδ(t) = 1 ∀ t. Second, we removed the proliferation term entirely by definingδ(t) = 0 ∀ t.

Cell to cell transmission
In addition to our original model including a free variable for viral load, we also explored the capacity of cell-to-cell transmission models to explain the available data. First, we formulated the standard mass action model as follows: for some saturation constant, s, andδ(t) is as defined in the main text. Second, we explored an alternative cell-to-cell transmission model that includes a correction factor, f c (I), to account for spread within a spatially confined region [5]. Such a model would be more applicable to dynamics occurring in a stationary environment, such as the lymphoid tissues, rather than that of the blood. The equations for this model are identical to Eqs S1-S5, except the transmission parameter β is multiplied by f c (I), where m is the number of neighbors of each cell (assuming spread within a two-dimensional lattice); and z, a, b are constants defined to maintain smoothness of f c ∀ I (see Kumberger et al. (2018) for further information [5]). For both these cell-to-cell transmission models, the viral load data were used to fit the infected cell compartment (I), the total lymphocyte data to fit the sum of the model lymphocyte compartments (S + I + A), and the MV-specific T cell data to fit the activated T cell compartment (A). We fit the model over a range of possible values for m (4, 6, and 8), and constant values for z, a, b were assigned following Kumberger et al.

Fitting procedure
For each individual and model, let n X be the number of measurements of variable X, where X is either total lymphocyte counts (L), viral load (V ), or MV-specific activated T cells (A). The residual sum of squares between the log-transformed data (D) and model predictions (M ) is then given by As outlined in Hogan et al. (2015), assuming normally-distributed residuals gives the following log-likelihood expression for X The maximum likelihood estimate of σ 2 X is then obtained by setting December 6, 2018 3/5 Substituting this into Eq S6 then gives and combining all three measurement variables gives the full log-likelihood

Uncertainty and sensitivity analyses
Uncertainty analysis was carried out using Latin Hypercube sampling (LHS), which assesses the variability in model outcome caused by changes in input parameters. The general idea is to randomly draw values for all input parameters from predefined probability distributions without replacement. Each sample of new draws (spanning the entire parameter set) is then used to re-simulate the model and calculate the new output, usually defined as some representative summary measure (e.g. peak viral load, total viral load, cumulative number of infection events etc). For N new parameter sets, one then has a distribution of N new outputs that can be used to assess model variability arising from parameter uncertainty. Previous work suggests N should satisfy N > 4/3K, where K is the number of input parameters [6].
In this study we used total viral load to describe model output, given its relevance to the experimental simulations comparing drivers of viral clearance. To explore model predictions regarding viral resurgence, we conducted additional analysis using the time to secondary viral rebound (measured in days post infection) as a supplemental measure. To avoid setting prior assumptions on parameter values, we defined uniform distributions for all input parameters; we also chose N = 100 (which satisfies 100 > 4/3 × 10 = 13.33). Note that the scaling factor, ψ, was not included in the analysis as it does not appear in the model equations and does not directly impact the resulting simulations. Bounds on fitted parameter distributions were defined as the maximum and minimum values of the corresponding fitted estimates across all individuals; bounds on fixed parameters were defined as follows: c ∈ (1, 6), d ∈ (1/80, 1/20), δ ∈ (1/4, 1). Widening these bounds did not impact our overall conclusions. All simulations were performed using the target cell and T cell model described by Eqs 1-6.
To assess model sensitivity to individual parameters, the above uncertainty analysis was extended using partial rank correlation (PRC) methods. In brief, the PRC coefficient of each input parameter was determined by correlating the LHS draws of that parameter with the corresponding model output (whilst controlling for all other parameters). Algorithms for calculating PRC coefficients are described in detail elsewhere (e.g. [6]). The sign of the PRC coefficient indicates the direction of the impact on model outcome (e.g. a positive correlation coefficient indicates an overall increase in the summary measure), whereas the magnitude indicates the strength of the impact. This technique therefore highlights the parameters with the greatest influence on model outcome. In summary, although the described LHS and PRC methods sample a range of possible model outcomes, rather than providing absolute bounds on the fitted model trajectories, they are powerful tools for assessing and visualizing the impact of parameter uncertainty on model predictions [6].