On the role of CD 8 + T cells in determining recovery time from influenza virus infection
Pengxing Cao, Zhongfang Wang, Ada W. C. Yan, Jodie McVernon, Jianqing Xu, Jane M. Heffernan, Katherine Kedzierska, James M. McCaw
aa r X i v : . [ q - b i o . CB ] S e p On the role of CD8 + T cells in determiningrecovery time from influenza virus infection
Pengxing Cao , Zhongfang Wang , Ada W. C. Yan , Jodie McVernon ,Jianqing Xu , Jane M. Heffernan , Katherine Kedzierska , and James M.McCaw ∗ School of Mathematics and Statistics, The University of Melbourne,Melbourne, Australia. Department of Microbiology and Immunology, University of Melbourne, atthe Peter Doherty Institute for Infection and Immunity, Parkville, Victoria,Australia Shanghai Public Health Clinical Center and Institutes of BiomedicalSciences, Key Laboratory of Medical Molecular Virology of Ministry ofEducation/Health, Shanghai Medical College, Fudan University, Shanghai,China Centre for Epidemiology and Biostatistics, Melbourne School of Populationand Global Health, The University of Melbourne, Melbourne, Australia. Modelling and Simulation, Infection and Immunity Theme, MurdochChildrens Research Institute, The Royal Children’s Hospital, Parkville,Victoria, Australia. Modelling Infection and Immunity Lab, Centre for Disease Modelling, YorkInstitute for Health Research, York University, Toronto, Ontario, Canada. ∗ Correspondence: [email protected] bstract
Myriad experiments have identified an important role for CD8 + T cell response mechanisms indetermining recovery from influenza A virus infection. Animal models of influenza infectionfurther implicate multiple elements of the immune response in defining the dynamical charac-teristics of viral infection. To date, influenza virus models, while capturing particular aspectsof the natural infection history, have been unable to reproduce the full gamut of observed viralkinetic behaviour in a single coherent framework. Here, we introduce a mathematical modelof influenza viral dynamics incorporating all major immune components (innate, humoral andcellular) and explore its properties with a particular emphasis on the role of cellular immunity.Calibrated against a range of murine data, our model is capable of recapitulating observed viralkinetics from a multitude of experiments. Importantly, the model predicts a robust exponentialrelationship between the level of effector CD8 + T cells and recovery time, whereby recoverytime rapidly decreases to a fixed minimum recovery time with an increasing level of effec-tor CD8 + T cells. We find support for this relationship in recent clinical data from influenzaA(H7N9) hospitalised patients. The exponential relationship implies that people with a lowerlevel of naive CD8 + T cells may receive significantly more benefit from induction of additionaleffector CD8 + T cells arising from immunological memory, itself established through eitherprevious viral infection or T cell-based vaccines.2 ntroduction
Invasion of influenza virus into a host’s upper respiratory tract leads to infection of healthyepithelial cells and subsequent production of progeny virions (1). Infection also triggers avariety of immune responses. In the early stage of infection a temporary non-specific response(innate immunity) contributes to the rapid control of viral growth while in the late stage ofinfection, the adaptive immune response dominates viral clearance (2). The early immuneresponse involves production of antiviral cytokines and cells, e.g. type 1 interferon (IFN) andnatural killer cells (NK cells), and is independent of virus type (3, 4, 5, 6, 7). In the special caseof a first infection in a naive host, the adaptive immune response, mediated by the differentiationof naive T cells and B cells and subsequent production of virus-specific T cells and antibodies(8, 2), leads to not only a prolonged killing of infected cells and virus but also the formationof memory cells which can generate a rapid immune response to secondary infection with thesame virus (9, 10).CD8 + T cells, which form a major component of adaptive immunity, play an important rolein efficient viral clearance (11). However, available evidence suggests they are unable to clearvirus in the absence of antibodies (12, 13) except in hosts with a very high level of pre-existingnaive or memory CD8 + T cells (14, 15, 16). Some studies indicate that depletion of CD8 + Tcells could decrease the viral clearance rate and thus prolong the duration of infection (17, 18,19, 20). Furthermore, a recent study of human A(H7N9) hospitalized patients has implicatedthe number of effector CD8 + T cells as an important driver of the duration of infection (21).This diverse experimental and clinical data, sourced from a number of host-species, indicatesthat timely activation and elevation of CD8 + T cell levels may play a major role in the rapidand successful clearance of influenza virus from the host. These observations motivate ourmodeling study of the role of CD8 + T cells in influenza virus clearance.Viral dynamics models have been extensively applied to the investigation of the antiviralmechanisms of CD8 + T cell immunity against a range of pathogens, with major contributionsfor chronic infections such as HIV/SIV (22, 23, 24, 25, 26, 27), HTLV-I (28) and chronicLCMV (29, 30). However, for acute infections such as measles (31) and influenza (32, 33, 34,35, 36, 37, 38, 39, 40, 41, 42), highly dynamical interactions between the viral load and theimmune response occur within a very short time window, presenting new challenges for thedevelopment of models incorporating CD8 + T cell immunity.Existing influenza viral dynamics models, introduced to study specific aspects of influenzainfection, are limited in their ability to capture all major aspects of the natural history of in-fection, hindering their use in studying the role of CD8 + T cells in viral clearance. Somemodels show a severe depletion of target cells (i.e healthy epithelial cells susceptible to viralinfection) after viral infection (34, 36, 37, 38, 40). Depletion may be due to either infection orimmune-mediated protection. Either way, these models are arguably incompatible with recentevidence that the host is susceptible to re-infection with a second strain of influenza a shortperiod following primary exposure (43). Furthermore, as reviewed by Dobrovolny et al. (39),target cell depletion in these models strongly limits viral expansion so that virus can be effec-tively controlled or cleared at early stage of infection even in the absence of adaptive immunity,which contradicts the experimental finding that influenza virus remains elevated in the absenceof adaptive immune response (44). While a few models do avoid target cell depletion (32, 33),they either assume immediate replenishment of target cells (32) or a slow rate of virus invasioninto target cells resulting in a much delayed peak of virus titer at day 5 post-infection (ratherthan the observed peak at day 2) (33). Moreover, models with missing or unspecified major im-mune components, e.g. no innate immunity (24, 25, 36, 38), no antibodies (24, 25, 33, 41, 42)3r unspecified adaptive immunity (40), also indicate the need for further model development.For an in-depth review of the current virus dynamics literature on influenza, we refer the readerto the excellent article by Dobrovolny et al. (39).
IFN (F)Targetcells (T)
Virus(V)
Antibody(A) p r o d u c e Infectedcells (I) p r o d u c e naiveB cells p r o d u c e naive CD8 + T cells stimulate proliferation &di (cid:1) erentiation s t i m u l a t e proliferation &di (cid:0) erentiation e (cid:2) ector CD8 + T cells plasmacells n e u t r a li s e k ill di (cid:3) erentantivirale (cid:4) ects Figure 1: Schematic diagram showing the major components of viral infection and the immuneresponse. Infection starts when virus binds to healthy epithelial cells (target cells). Infectedcells release new virus and produce cytokines such as IFN. IFN is a major driver of innateimmunity, responsible for effective control of rapid viral growth and expansion. Virus furtherstimulates naive CD8 + T cells and B cells to produce effector CD8 + T cells and antibodies,responsible for final clearance of virus.In this paper, we construct a within-host model of influenza viral dynamics in naive (i.e.previously unexposed) hosts that incorporates the major components of both innate and adap-tive immunity and use it to investigate the role of CD8 + T cells in influenza viral clearance.The model is calibrated against a set of published murine data from Miao et al. (38) and is thenvalidated through demonstration of its ability to qualitatively reproduce a range of publisheddata from immune-knockout experiments (12, 17, 18, 44, 13, 38). Using the model, we find thatthe recovery time — defined to be the time when virus titer first drops below a chosen thresholdin the (deterministic) model — is negatively correlated with the level of effector CD8 + T cellsin an approximately exponential manner. To the best of our knowledge, this relationship, withsupport in both H3N2-infected mice and H7N9-infected humans (21), has not been previouslyidentified. The exponential relationship between CD8 + T cell level and recovery time is shownto be remarkably robust to variation in a number of key parameters, such as viral productionrate, IFN production rate, delay of effector CD8 + T cell production and the level of antibodies.Moreover, using the model, we predict that people with a lower level of naive CD8 + T cellsmay receive significantly more benefit from induction of additional effector CD8 + T cells. Suchproduction, arising from immunological memory, may be established through either previous4iral infection or T cell-based vaccines.
Methods
The model
The model of primary viral infection is a coupled system of ordinary and delay differentialequations, consisting of three major components (see Fig. 1 for a schematic diagram). Eqs. 1–3describe the process of infection of target cells by influenza virus and are a major component inalmost all models of virus dynamics in the literature. Eqs. 4 and 5 model IFN-mediated innateimmunity (45, 46). Thirdly, adaptive immunity including CD8 + T cells and B cell-producedantibodies for killing infected cells and neutralizing influenza virus respectively are describedby Eqs. 6–11. dVdt = p V I − d V V − k S VA S − k L VA L − b V T , (1) dTdt = g T ( T + R )( − T + R + IT ) − b ′ V T + r R − f F T , (2) dIdt = b ′ V T − d I I − k N IF − k E IE , (3) dFdt = p F I − d F F , (4) dRdt = f FT − r R , (5) dC n dt = − b Cn ( VV + h C ) C n , (6) dEdt = b Cn ( V ( t − t C ) V ( t − t C ) + h C ) C n ( t − t C ) e ( p C t C ) − d E E , (7) dB n dt = − b Bn ( VV + h B ) B n , (8) dPdt = b Bn ( V ( t − t B ) V ( t − t B ) + h B ) B n ( t − t B ) e ( p B t B ) − d P P , (9) dA S dt = p S P − d S A S , (10) dA L dt = p L P − d L A L . (11)In further detail, Eq. 1 indicates that the change in viral load ( V ) is controlled by fourfactors: the production term ( p V I ) in which virions are produced by infected cells ( I ) at a rate p V (47, 37, 45); the viral natural decay/clearance ( d V V ) with a decay rate of d V ; the viralneutralisation terms ( k S VA S and k L VA L ) by antibodies (both a short-lived antibody response A S driven by, e.g. IgM, and a longer-lived antibody response A L driven by, e.g. IgG and IgA(12, 38)), and a consumption term ( b V T ) due to binding to and infection of target cells ( T ).In Eq. 2, the term g T ( T + R )( − ( T + R + I ) / T ) models logistic regrowth of the target cellpool (46). Both target cells ( T ) and resistant cells ( R , those protected due to IFN-inducedantiviral effect) can produce new target cells, with a net growth rate proportional to the severityof infection, 1 − ( T + R + I ) / T (i.e. the fraction of dead cells). T is the initial number of target5ells and the maximum value for the target cell pool (34). Target cells ( T ) are consumed byvirus ( V ) due to binding ( b ′ V T ), the same process as b V T . Note that b and b ′ have differentmeasurement units due to different units for viral load ( V ) and infected cells ( I ). As alreadymentioned, the innate response may trigger target cells ( T ) to become resistant ( R ) to virus,at rate f F T . Resistant cells lose protection at a rate r (45). This process also governs theevolution of virus-resistant cells ( R ) in Eq. 5.Eq. 3 describes the change of infected cells ( I ). They increase due to the infection of targetcells by virus ( b ′ V T ) and die at a (basal) rate d I . Two components of the immune responseincrease the rate of killing of infected cells. IFN-activated NK cells kill infected cells at a rate k N IF (48, 6, 45, 46). Effector CD8 + T cells ( E ) — produced through differentiation from naiveCD8 + T cells C n in Eq. 6 — kill at a rate k E IE . Of note our previous work has demonstratedthat models of the innate response containing only IFN-induced resistance for target cells (state R ; Eq. 5), while able to maintain a population of healthy uninfected cells, still control viralkinetics through target cell depletion, and therefore cannot reproduce viral re-exposure data(46, 43). Given our interest in analysing a model that prevents target cell depletion, inclusionof IFN-activated NK cells (term k N IF ) is an essential part of the model construction.Eq. 4 models the innate response, as mediated by IFN ( F ). IFN is produced by infectedcells at a rate p F and decays at a rate d F (46).Eq. 6 models stimulation of naive CD8 + T cells ( C n ) into the proliferation/differentiationprocess by virus at a rate b Cn V / ( V + h C ) ), where b Cn is the maximum stimulation rate and h C indicates the viral load ( V ) at which half of the stimulation rate is achieved. Note that thisformulation does not capture the process of antigen presentation and CD8 + T cell activation,but rather is a simple way to establish the essential coupling between the viral load and the rateof CD8 + T cell activation in the model (49). In Eq. 7, the production of effector CD8 + T cells( E ) is assumed to be an “advection flux” induced by a delayed virus-stimulation of naive CD8 + T cells (the first term on the righthand side of Eq. 7). The delayed variables, V ( t − t C ) and C n ( t − t C ) , equal zero when t < t C . The introduction of the delay t C is to phenomenologicallymodel the delay induced by both naive CD8 + T cell proliferation/differentiation and effectorCD8 + T cell migration and localization to the site of infection for antiviral action (50, 51,42). The delay also captures the experimental finding that naive CD8 + T cells continue todifferentiate into effector T cells in the absence of ongoing antigenic stimulation (49, 52). Themultiplication factor e p C t C indicates the number of effector CD8 + T cells produced from onenaive CD8 + T cell, where p C is the average effector CD8 + T cell production rate over thedelay period t C . The exponential form of the multiplication factor is derived based on theassumption that cell differentiation and proliferation follows a first-order advection–reactionequation. Effector CD8 + T cells decay at a rate d E .Similar to CD8 + T cells, Eqs. 8 and 9 model the proliferation/differentiation of naive Bcells, stimulated by virus presentation at rate b Bn V / ( V + h B ) . Stimulation subsequently leadsto production of plasma B cells ( P ) after a delay t B . The multiplication factor e p B t B indicatesthe number of plasma B cells produced from one naive B cell, where p B is the production rate.Plasma B cells secrete antibodies, which exhibit two types of profiles in terms of experimentalobservation: a short-lived profile (e.g. IgM lasting from about day 5 to day 20 post-infection)and a longer-lived profile (e.g. IgG and IgA lasting weeks to months) (12, 38). These twoantibody responses are modeled by Eqs. 10 and 11 wherein different rates of production ( P S and P L ) and consumption ( d S and d L ) are assumed.6 odel parameters and simulation The model contains 11 equations and 30 parameters (see Table 1). This represents a seriouschallenge in terms of parameter estimation, and clearly prevents a straightforward applicationof standard statistical techniques. To reduce uncertainty, a number of parameters were takendirectly from the literature, as per the citations in Table. 1. The rest were estimated (as indi-cated in Table 1) by calibrating the model against the published data from Miao et al. (38)who measured viral titer, CD8 + T cell counts and IgM and IgG antibodies in laboratory mice(exhibiting a full immune response) over time during primary influenza H3N2 virus infection(see (38) for a detailed description of the experiment). The approach to estimating the param-eters based on Miao et al. ’s data is provided in the
Supplementary Material and the estimatedparameter values are given in Table 1. Note that the data were presented in scatter plots in theoriginal paper (38), while we presented the data here in Mean ± SD at each data collectiontime point for a direct comparison with our mean-field mathematical model.For model simulation, the initial condition is set to be ( V , T , I , F , R , C n , E , B n , P , A S , A L ) =( V , T , , , , , , , , , ) unless otherwise specified. The initial target cell number ( T )was estimated by Petrie et al. (53). We estimate that of order 100 cells (resident in the spleen)are able to respond to viral infection ( C n ) (personal communication, N. LaGruta, Monash Uni-versity, Australia). Note that 100 naive CD8 + T cells might underestimate the actual number ofnaive precursors that could respond to all the epitopes contained within the virus but does notqualitatively alter the model dynamics and predictions (see
Results where the naive CD8 + Tcell number is varied between 0 to 200). In the absence of further data, we also use this value forthe initial naive B cell number ( B n ), but again this choice does not qualitatively alter the modelpredictions. The numerical method and code (implemented in MATLAB, version R2014b, theMathWorks, Natick, MA) for solving the model are provided in the Supplementary Material . Analysis of clinical influenza A(H7N9) data
Clinical influenza A(H7N9) patient data was used to test our model predictions on the relation-ship between CD8 + T cell number and recovery time. The data was collected from 12 survivingpatients infected with H7N9 virus during the first wave of infection in China in 2013 (raw datais provided in
Dataset S1 ; see the paper of Wang et al. (21) for details of data collection; thisstudy was reviewed and approved by the SHAPHC Ethics Committee). Note that the clinicaldata were scarce for some patients. For those patients, we have assumed that the available dataare representative of the unobserved values in the neighboring time period. For each patient,we took the average IFN g + CD8 + T cell number in 10 peripheral blood mononuclear cells(PBMC) for the period from day 8 to day 22 (or the recovery day if it comes earlier) post-admission as a measure of the effector CD8 + T cell level. This period was chosen a priori as itroughly matches the duration of the CD8 + T cell profile and clinical samples were frequentlycollected in this period. The average CD8 + T cell count was given by the ratio of the total areaunder the data points (using trapezoidal integration) to the number of days from day 8 to day 22(or the recovery day if it comes earlier). For those patients for which samples at days 8 and/or22 were missing we specified the average CD8 + T cell level at the missing time point to beequal to the value from the nearest sampled time available.7
10 20 l o g ( V ) ( E I D / m l ) -202468 0 10 20 30 C D + T c e ll s ×10 I g M ( pg / m l ) time (days) I g G ( pg / m l ) Figure 2: The model with estimated parameters (solid curves) captures the murine data fromthe paper by Miao et al. (38). The data is shown with error bars (Mean ± SD). Note that due tothe limit of detection for the viral load (occurring after 10 days post-infection as seen in viralload data), the last three data points in the upper-left panel were not taken into considerationfor model fitting.
Results
Model properties and reproduction of published experimental data
We first analyze the model behavior in order to establish a clear understanding of the modeldynamics. Fig. 2 shows solutions (time-series) for the model compartments (viral load, CD8 + T cells and IgM and IgG antibody) calibrated against the murine data from (38). Solutionsfor the remaining model compartments are shown in Fig. 3. The model (with both innate andadaptive components active) prevents the depletion of target cells (see Fig. 3 wherein over 50%of target cells remain during infection) and results in a minor loss of just 10–20% of healthyepithelial cells (i.e. the sum of target cells ( T ) and virus-resistant cells ( R ); see SupplementaryFig. S1). The primary driver for the maintenance of the target cell pool during acute viralinfection is a timely activation of the innate immune response, and in particular the natural killercells (Supplementary Fig. S2). Since we have previously shown that the target-cell limitedmodel (even with the resistant cell compartment) is unable to reproduce observations fromheterologous re-exposure experiments (43, 46), our model improves upon previous models8
10 20 T ( c e ll s ) × I ( c e ll s ) × R ( c e ll s ) × F C n ( c e ll s ) B n ( c e ll s ) time (days) P ( c e ll s ) Figure 3: Model solution with parameters given in Table 1. Time courses of viral load ( V ),effector CD8 + T cells ( E ), short-lived antibody response ( A S ) and long-lived antibody response( A L ) have been shown in Fig. 2.where viral clearance was only achieved through depletion of target cells (a typical solutionshown in Supplementary Fig. S2b). Importantly, our result is distinguished from that of Saenzet. al. (37), wherein the healthy cell population was similarly maintained, but primarily throughinduction of the virus-resistant state, thereby rendering that model incapable of capturing re-9nfection behavior as established in animal models (43).The modeled viral dynamics exhibits three phases, each dominated by the involvement ofdifferent elements of the immune responses (Fig. 4). Immediately following infection (0–2 dayspost-infection) and prior to the activation of the innate (and adaptive) immune responses, virusundergoes a rapid exponential growth (Fig. 4a). In the second phase (2–5 days post-infection),the innate immune response successfully limits viral growth (Fig. 4a). In the third phase (4–6days post-infection), adaptive immunity (antibodies and CD8 + T cells) is activated and viralload decreases rapidly, achieving clearance. Figs. 4b and 4c demonstrate the dominance ofthe different immune mechanisms at different phases. In Fig. 4b models with and withoutimmunity are indistinguishable until day 2 (shaded region), before diverging dramatically whenthe innate and then adaptive immune responses influence the dynamics. In Fig. 4c, modelswith and without an adaptive response only diverge at around day 4 as the adaptive responsebecomes active. We have further shown that this three-phase property is a robust feature of themodel, emergent from its mathematical structure and not a property of fine tuning of parameters(see Supplementary Fig. S3). Importantly, it clearly dissects the periods and effect of innateimmunity, extending on previous studies of viral infection phases where the innate immuneresponse was either ambiguous or ignored (12, 54, 38).As reviewed by Dobrovolny et al. (39), a number of in vivo studies have been performedto dissect the contributions of CD8 + T cells and antibodies (12, 17, 18, 44, 55). We use thefindings of these studies to validate our model, by testing how well it is able to reproduce theexperimental findings (without any further adjustment to parameters). Although the determi-nation of the role of CD8 + T cells is often hindered by co-inhibition of both CD8 + T cellsand the long-lived antibody response (e.g. using nude mice), it is consistently observed thatantibodies play a dominant role in final viral clearance while CD8 + T cells are primarily re-sponsible for the timely killing of infected cells and so indirectly contribute to an increased rateof removal of free virus towards the end of infection (17, 18, 13). Furthermore, experimentaldata demonstrate that a long-lived antibody response is crucial for achieving complete viralclearance, while short-lived antibodies are only capable of driving a transient decrease in viralload (12, 44). We find that our model (with parameters calibrated against Miao et al. ’s data(38)) is able to reproduce these observations: • virus can rebound in the absence of long-lived antibody response (see Fig. 5 and Supple-mentary Fig. S4). • both the CD8 + T cell response and short-lived antibody response only facilitate a fasterviral clearance, and are incapable of achieving clearance in the absence of long-livedantibody response (see Fig. 5 and Supplementary Fig. S4). • a lower level of CD8 + T cells (modulated by a decreased level of initial naive CD8 + Tcells, C n ) significantly prolongs the viral clearance (see Supplementary Fig. S4).In addition, the model also predicts a rapid depletion of naive CD8 + T cells after primaryinfection (see Fig. 3), which represents a full recruitment of naive CD8 + T cell precursors.This result may be associated with the experimental evidence suggesting a strong correlationbetween the naive CD8 + T cell precursor frequencies and effector CD8 + T cell magnitudesfor different pMHC-specific T cell populations (56). Note that in Fig. 5 no adjustments tothe model (e.g. to the vertical scale) were made; its behavior is completely determined by thecalibration to the aforementioned murine data (38) and so these findings represent a (successful)prediction of the model. 10 ime (days)0 10 20 l o g ( V ) ( E I D / m l ) full immuneno immunity time (days) l o g ( V ) ( E I D / m l ) full immuneIFN only (days) l o g ( v i r a l l o a d ) e xponential growth innate immune response adaptive immune response (a)(b) (c) Figure 4: The model solution exhibits three-phase behavior following influenza virus infection.Panel (a) schematically represents these three-phases of behavior based on involvement of im-mune responses. Following infection, virus first undergoes rapid exponential growth before theinnate immune response is activated (on day 2). Innate immunity controls rapid viral expansionto form a plateau in viral load. Adaptive immunity starts to take effect 4–6 days post-infectionand is responsible for viral clearance. Panels (b) and (c) demonstrate that the model dynam-ics follow the three-phase theory. Viral kinetics in the presence of a full immune response isshown by the solid line (in (b) and (c)). In panel (b), the dashed line shows viral kinetics inthe complete absence of immunity (innate and adaptive; by letting p F = b Cn = b Bn = b Cn = b Bn =
0; innate immunity remains active). The trajectories overlap prior tothe activation of the adaptive response. The shaded region highlights the second phase (innateresponse). Note: changes in model parameters shifts where the three phases occur, but doesnot alter the underlying three-phase structure. i.e. existence of the three phases is robust tovariation in parameters (see
Supplementary Material and Supplementary Fig. S3 in particular).In summary, we have demonstrated that our model—with parameters calibrated againstmurine data (38)—exhibits three important phases characterized by the involvement of variousimmune responses. Advancing on previous models, our model does not rely on target celldepletion, and successfully reproduces a multitude of behavior from knockout experimentswhere particular components of the adaptive immune response were removed. This providesus with some confidence that each of the major components of the immune response has beencaptured adequately by our model, allowing us to now make predictions on the effect of thecellular adaptive response on viral clearance. 11 ime (days) p u l m o n a r yv i r u s t i t e r ( P F U ) l o g f ull immune I g A-(cid:5) g G(cid:6)(cid:7) g M(cid:8), (cid:9) g (cid:10)(cid:11) time (days) (cid:12) V i r u s l o g ( E(cid:13)D ) (cid:14) (cid:15) normal B(cid:16)L(cid:17)/c mi (cid:18) e ( (cid:19) ull immune)nude mi (cid:20) e ( C(cid:21) T(cid:22)(cid:23) (cid:24) g (cid:25)(cid:26)(cid:27) (cid:28) g (cid:29)(cid:30) )nude mi (cid:31) e + external b s time (days) ! " l o g ( )( $%& ’ m l ) ( ) * ull immune .0 g g =>?@ F HJK N g OPQ R g SUW external X g Y time (days) Z l o g ( [ )( \]^ _ m l ) ‘ a d ull immuneno long e lived anti g ody response (a)(b) Figure 5: Consistency between mice data (left panels) and model results (right panels) showsthat short-lived antibody response (e.g. IgM) is only capable of generating a transient decreasein viral load while long-lived antibody response (e.g. IgA) plays a more dominant role inlate-phase viral clearance. (a) Data is from the paper of Kris et al. (44). Normal or nudeBALB/c mice were infected with H3N2 virus. External antibodies were given at day 5 andhad waned by about day 21. The model simulation mimics the passive antibody input byintroducing an extra amount of IgM (in addition to antigen-stimulated IgM), whose time coursefaithfully follows the experimental measurement (see Fig. 2a in the paper of Kris et al. (44)).“CD8 + T-, IgA-, IgG-” was modeled by letting b Cn = p L =
0. “External IgM” (inaddition to the IgM produced by plasma cells) was modeled by adding a new term − k S VA e ( t ) to Eq. 1 where A e ( t ) follows a piecewise function A e ( t ) = t < A e ( t ) = ( t − ) for t ∈ [ , ) , A e ( t ) = − ( t − ) for t ∈ [ , ) and A e ( t ) = t ≥
21. (b) Data is fromthe paper of Iwasaki et al. (12). The data indicates that the long-lasting IgA response, butnot the long-lasting IgG response or the short-lasting IgM response, is necessary for successfulviral clearance. “No long-lived antibody response” was modeled by letting p L =
0. Note thatMiao et al. only measured IgM and IgG, but not IgA. As such, our model’s long-lived antibodyresponse was calibrated against IgG kinetics (see Fig. 2). Therefore, we emphasise that we canonly investigate the relative contributions of short-lived and long-lived antibodies.
Dependence of the recovery time on the level of effector
CD8 + T cells
Having established that our model is (from a structural point of view) biologically plausibleand that our parameterization is capable of reproducing varied experimental data under differ-ent immune conditions (i.e. knockout experiments), we now study how the cellular adaptiveresponse influences viral kinetics in detail. We focus on the key clinical outcome of recov-ery time, defined in the model as the time when viral titer first falls below 1 EID / ml, the12inimum value detected in relevant experiments (e.g. Fig. 2). naive CD8 + T cell number (cells) a v e r a g e C D + T c e ll s ( × c e ll s ) h i jkl m model simulationlinear n it naive op qr s ell num t er ( u ells)0 50 100 150 v r e w o v e r y t i m e ( d a y s ) x y z model simulationexponential { it average |} ~(cid:127) (cid:128) ells ( (cid:215) (cid:129) ells)0 1 (cid:130) r e (cid:131) o v e r y t i m e ( d a y s ) (cid:132) (cid:133) model simulation (varying (cid:134) n ) exponential (cid:135) it average (cid:136)(cid:137) (cid:138)(cid:139) (cid:140) ells ( (cid:141) (cid:142) ells)0 1 (cid:143) r e (cid:144) o v e r y t i m e ( d a y s ) (cid:145) (cid:146) model simulation (varying (cid:147) n )model simulation (varying (cid:148)(cid:149) ) (a) (d)(b)(c) Figure 6: The level of effector CD8 + T cells plays an important role in determining recoverytime. Recovery time is defined to be the time when viral load falls to 1 EID / ml. Panel (a)shows that the average effector CD8 + T cell number over days 6–20 is linearly related to thenaive CD8 + T cell number (i.e. C n ( ) ). Panel (b) shows that the recovery time is approximatelyexponentially related to the initial naive CD8 + T cell number. Combined, these results givepanel (c) wherein an approximately exponential relationship is observed between the averageCD8 + T cell number and recovery time, both of which are experimentally measurable. Notethat the exponential/linear fits shown in the figures are not generated by the viral dynamicsmodel but are used to indicate the trends (evident visually) in the model’s behavior. Panel (d)shows that varying the delay t C (in a similar way to that shown in Fig. S5 in the SupplementaryMaterial ), rather than the naive CD8 + T cell number, does not alter the exponential relationship.In panel (d), the crosses represent the results of varying t C and the empty circles are the sameas those in panel (c) for comparison.Time series of the viral load show that the recovery time decreases as the initial naive CD8 + T cell number ( C n ) increases (Supplementary Fig. S4). With that in mind, we now examinehow recovery time is associated with the clinically relevant measure of effector CD8 + T celllevel during viral infection. With an increasing initial level of naive CD8 + T cells, the averagelevel of effector CD8 + T cells over days 6–20 increases linearly (Fig. 6a), while the recoverytime decreases in an approximately exponential manner (Fig. 6b). Combining these two effectsgives rise to an approximately exponential relation between the level of effector CD8 + T cellsand recovery time (Fig. 6c). Note that the exponential/linear fits shown in the figures are simplyto aid in interpretation of the results. They are not generated by the viral dynamics model.13f varying the delay for naive CD8 + T cell activation and differentiation, t C , while keepingthe naive CD8 + T cell number fixed (at the default value of 100), we find that the average levelof effector CD8 + T cells is exponentially related to the delay, while the recovery time is depen-dent on the delay in a piecewise linear manner (see Supplementary Fig. S5). Nevertheless, thecombination still leads to an approximately exponential relationship between the level of effec-tor CD8 + T cells and recovery time (Supplementary Fig. S6), which is almost identical to thatof varying naive CD8 + T cells (Fig. 6d). We also examine the sensitivity of the exponentialrelationship to other model parameters generally accepted to be important in influencing themajor components of the system, such as the viral production rate p V , IFN production rate p F and naive B cell number. We find that the exponential relationship is robust to significant vari-ation in all of these parameters (see Supplementary Figs. S6 and S7 and Fig. 9). These resultssuggest that a higher level of effector CD8 + T cells is critical for early recovery, consistent withexperimental findings (57).
No. of CD8 + T cells in 10 PBMC ( × )0 1 2 3 4 r e c o v e r y t i m e ( da ys ) data (all excluding patient a79)data of patient a79exponential fit Figure 7: Clinical A(H7N9) data exhibits an exponential relationship between the averageCD8 + T cell number and recovery time. The x-axis is the average level of functional ef-fector CD8 + T cells (i.e. IFN g + CD8 + T cells) over the period from day 8 to day 22(or the recovery day if it comes earlier). Spearman’s rank correlation test indicates a sig-nificant negative correlation between the average CD8 + T cell numbers and recovery time( r = − . , p = . Dataset S1 ; discussedin the
Discussion ), all other data points (solid dots) are fitted by an offset exponential function y = . e − . x + . + T cell response is approximately 17.5356.Finally, and perhaps surprisingly given our model has been calibrated purely on data fromthe mouse, a strikingly similar relationship as shown in Fig. 6c is found in clinical data frominfluenza A(H7N9) virus-infected patients (Fig. 7). Excluding one patient (No. a79 in
DatasetS1 ; the exclusion is considered further in the
Discussion ), average IFN g + CD8 + T cells and re-covery time are negatively correlated (Spearman’s r = − . p = . + T cells and a strong saturation for relatively highCD8 + T levels, implying that even with a very high level of naive CD8 + T cells, recovery timecan not be reduced below a certain value (in this case, estimated to be approximately 17 days).Of course, the exponential relationship (i.e. the scale of CD8 + T cell level or recovery time),is only a qualitative one, as we have no way to determine the scaling between different x-axismeasurement units, nor adjust for particular host and/or viral factors that differ between the twoexperiments (i.e. H3N2-infection in the mouse (38) versus H7N9-infeciton in humans (21)).
Dependence of the recovery time on the level of memory
CD8 + T cells
In addition to naive CD8 + T cells, memory CD8 + T cells (established through previous viralinfection) may also significantly affect recovery time due to both their rapid activation uponantigen stimulus and faster replication rate (58, 59, 60, 61). To study the role of memory CD8 + T cells, we must first extend our model. As we are only concerned with how the presence ofmemory CD8 + T cells influences the dynamics, as opposed to the development of the memoryresponse itself, the model is modified in a straightforward manner through addition of twoadditional equations which describe memory CD8 + T cell ( C m ) proliferation/differentiation: dC m dt = − b Cm ( VV + h Cm ) C m , (12) dE m dt = b Cm ( V ( t − t Cm ) V ( t − t Cm ) + h Cm ) C m ( t − t Cm ) e ( p Cm t Cm ) − d E E m . (13)Accordingly the term k E IE in Eq. 3 is modified to k E I ( E + E m ) . The full model and detailson the choice of the additional parameters are provided in the Supplementary Material . Notethat the model component, C m , may include different populations of memory CD8 + T cells,including those directly specific to the virus and those stimulated by a different virus but whichprovide cross-protection (62, 63).Fig. 8a shows how the pre-existing memory CD8 + T cell number ( C m ) changes the ex-ponential relationship between naive CD8 + T cells and recovery time. Importantly, as thenumber of memory CD8 + T cells increases, the recovery time decreases for any level of naiveCD8 + T cells and the exponential relationship remains. The extent of reduction in the recoverytime for a relatively low level of naive CD8 + T cells is greater than that for a relatively highlevel of naive CD8 + T cells. This suggests that people with a lower level of naive CD8 + Tcells may benefit more through induction of memory CD8 + T cells, emphasizing the potentialimportance for taking prior population immunity into consideration when designing CD8 + Tcell-based vaccines (64).The above result is based on the assumption that the initial memory CD8 + T cell numberupon re-infection is independent of the number of naive CD8 + T cells available during theprevious infection. However, it has also been found that the stationary level of memory CD8 + T cells is usually maintained at about 5–10% of the maximum antigen-specific CD8 + T cellnumber during primary viral infection (8, 65). This indicates that people with a low naive CD8 + T cell number may also develop a low level of memory CD8 + T cells following infection. Inconsequence, such individuals may be relatively more susceptible to viral re-infection (66).This alternative and arguably more realistic relationship between the numbers of naive andmemory CD8 + T cells is simulated in Fig. 8b where memory CD8 + T cell levels are set to 5%of the maximum of the effector CD8 + T cell level. Results suggest that, upon viral re-infection,pre-existing memory CD8 + T cells are able to significantly improve recovery time except for15 aive CD8 + T cell number (cells) r e c o v e r y t i m e ( d a y s ) C m (cid:150) (cid:151) m (cid:152) (cid:153) m (cid:154) (cid:155) (cid:156) m (cid:157) (cid:158) m (cid:159) naive (cid:160)¡ ¢£ ⁄ ell num ¥ er ( ƒ ells)0 10 § r e ¤ o v e r y t i m e ( d a y s ) ' “ « no memory ‹› fifl (cid:176) ells with memory (1 – )with memory (5 † ) (a) (b) Figure 8: Effects of naive and memory CD8 + T cells on viral clearance. Recovery time isdefined to be the time when the viral load falls to 1 EID / ml. Panel (a) demonstrates thatvarying the number of memory CD8 + T cells ( C m ) reduces the recovery time for any naiveCD8 + T cell number (i.e. C n ( ) ). Note that saturation is observed for C m >
100 where therecovery time is about 6 days, independent of the naive cell numbers. Panel (b) demonstrateshow the presence of pre-existing memory CD8 + T cells (solid dots) leads to a shorter recoverytime when compared to the case where no memory CD8 + T cells are established (open circles).Note the time scale difference in panels (a) and (b). This simulation is based on the assumptionthat the level of pre-existing memory CD8 + T cells is assumed to be either 1% or 5% (asindicated in the legend) of the maximum effector CD8 + T cell number due to primary viralinfection. The memory cell number (which is not shown in this figure) is about 30 time asmany as the naive cell number shown in the figure, i.e. 30 naive cells result in about 900memory cells before re-infection.in hosts with a very low level of naive CD8 + T cells (Fig. 8b). This is in accordance withthe assumption that a smaller naive pool leads to a smaller memory pool and in turn a weakershortening in recovery time. Although the model suggests that the failure of memory CD8 + T cells to protect the host is unlikely to be observed (because of the approximately 30 foldincrease in the size of the memory pool relative to the naive pool), the failure range may beincreased if the memory pool size is much smaller (modulated by, say, changing 5% to 1% inthe model). Therefore, for people with a low naive CD8 + T cell number, the level of memoryCD8 + T cells may be insufficient and prior infection may provide very limited benefit, furtheremphasizing the opportunity for novel vaccines that are able to induce a strong memory CD8 + T cell response to improve clinical outcomes.
Dependence of the recovery time on antibody level
Antibodies appear at a similar time as effector CD8 + T cells during influenza viral infectionand may enhance the reduction in the recovery time in addition to CD8 + T cells. By varyingthe naive B cell number B n (as a convenient, but by no means unique, way to influence antibodylevel), we find that increasing the antibody level shortens the recovery time regardless of theinitial naive CD8 + T cell number, leaving the exponential relation largely intact (Fig. 9). Aslight saturation occurs for the case in which levels of both naive B cells and CD8 + T cells arelow. Moreover, variation in naive B cell number also results in a wider variation in recoverytime for a lower naive CD8 + T cell level, suggesting that people with a lower level of naive16 aive CD8 + T cell number (cells)0 50 100 150 200 r e c o v e r y t i m e ( da ys ) B n = 60B n = 90B n = 120B n = 150 Figure 9: Influence of antibody level on the relationship between the naive CD8 + T cell numberand recovery time. Recovery time is defined to be the time when viral load falls to 1 EID / ml.Different antibody levels are simulated by varying the initial number of naive B cells (i.e. B n at t = + T cells may, once again, receive a more significant benefit (in terms of recovery time)through effective induction of an antibody response via vaccination.
Discussion
In this paper, we have studied the role of CD8 + T cells in clearing influenza virus from thehost using a viral dynamics model. The model was calibrated on a set of published murine datafrom Miao et al. (38) and has been further shown to be able to reproduce a range of publisheddata from experiments with different immune components knocked out. By avoiding target-cell depletion, our model is also compatible with re-infection data (43), providing a strongplatform on which to examine the role of CD8 + T cells in determining recovery time frominfection. Our primary finding is that the time of recovery from influenza viral infection isnegatively correlated with the level of effector CD8 + T cells in an approximately exponentialmanner. This robust property of infection has been identified from the model when calibratedagainst influenza A(H3N2) infection data in mice (38), but also observed in clinical case seriesof influenza A(H7N9) infection in humans (Fig. 7) (21). Our findings, in conjunction withconclusions on the potential role for a T cell vaccine that stimulates and/or boosts the memoryresponse, suggest new directions for research in both non-human species and further studies inhumans on the association between CD8 + T cell levels and clinical outcomes. Further research,including detailed statistical fitting of our model to an extensive panel of infection data (as yetunavailable) from human and non-human species, is required to establish the generality of theserelationships and provide quantitative insights for specific viruses in relevant hosts.The non-linear relationship between effector CD8 + T cell level and recovery time may beuseful in clinical treatment. The saturated property of the relation implies that a linear increase17n the effector CD8 + T cell level may result in diminishing incremental improvements in patientrecovery times. With evidence of a possible age-dependent loss of naive T cells (67, 68, 69),our model results imply that boosting the CD8 + T cell response via T cell vaccination maybe particularly useful for those with insufficient naive CD8 + T cells. The population-levelconsequences of such boosting strategies, while beyond the scope of this work, have previouslybeen considered by the authors (64).We also investigated the effect of memory CD8 + T cell level on viral clearance and foundunsurprisingly that a high pre-existing level of memory CD8 + T cells was always beneficial.However, our results suggest that pre-existing memory CD8 + T cells may be particularly ben-eficial for certain groups of people. For example, if the memory CD8 + T cell number inducedby viral infection or vaccination is assumed to be relatively constant for everyone, people withless naive CD8 + T cells would benefit more upon viral re-infection (see Fig. 8a). On the otherhand, if assuming pre-existing memory CD8 + T cell number is positively correlated with thenumber of naive CD8 + T cells (simulated in Fig. 8b), people with more naive CD8 + T cellswould benefit more upon viral re-infection. Emerging evidence suggests that the relationshipbetween the level of memory CD8 + T cells and naive precursor frequencies is likely to bedeeply complicated (56, 70, 71, 72). In that context, our model predictions emphasize the im-portance for further research in this area, and the necessity to take prior population immunityinto consideration when designing CD8 + T cell vaccines (64).We modeled both short-lived and long-lived antibody responses. Experimental data andmodel predictions consistently show that the short-lived antibody response results in a tem-porary reduction in virus level whereas the long-lived antibody response is responsible forcomplete viral clearance (Fig. 5). We emphasize here that although the model is able to capturethe observed short-lived and long-lived antibody responses (in order to study the virus-immuneresponse interactions), it is not designed to investigate the mechanisms inducing different an-tibody responses. The observed difference in antibody decay profile may be a result of manyfactors including the life times of different antibody-secreting cell types (73, 74), different an-tibody life times (75) and antibody consumption through neutralizing free virions. Detailedstudy of these phenomena requires a more detailed model and associated data for parameterestimation and model validation, and is thus left for future work. Similarly, CD4 + T cells arealso known to perform a variety of functions in the development of immunity, such as facil-itation of the formation and optimal recall of CD8 + T cells or even direct killing of infectedcells during viral infection (9, 10, 76, 77). Their depletion due to, say, HIV infection has alsobeen associated with more severe clinical outcomes following influenza infection (78). Someof the major functions of CD4 + T cells may be considered to be implicitly modeled throughrelevant parameters such as the rate of recall of memory CD8 + T cells (modeled by the delay t Cm ) in our extended model which includes memory CD8 + T cells. However, a detailed viraldynamics study of the role of CD4 + T cells in influenza infection, including in HIV infectedpatients with depleted CD4 + T cells, remains an open and important challenge.In a recent theoretical study, it was found that spatial heterogeneity in the T cell distributionmay influence viral clearance (42). Resident CD8 + T cells in the lungs have a more directand significant effect on timely viral clearance than do naive and memory pools resident inlymph nodes. Although this factor has been partially taken into consideration in our model byintroducing a delay for naive/memory CD8 + T cells, lack of explicit modeling of the spatialdynamics limits a direct application of our model to investigate these spatial effects.Finally, as noted in the results, one of the influenza A(H7N9)-infected patients (patienta79) was not included in our analysis of the clinical data (Fig. 7). Although our model suggestssome possibilities for the source of variation due to possible variation in parameter values, large18ariations in recovery time are only expected to occur for relatively low levels of naive CD8 + Tcells, nominally incompatible with this patient’s moderate CD8 + T cell response but a relativelylong recovery time. However, we note that IFN g + CD8 + T cell counts for this patient wereonly collected at days 10 and 23, and that the count at day 10 was particularly low and muchlower than that at day 23 (see
Dataset S1 ). We suspect that delayed, rather than weakened,production (to at least day 10) of the IFN g + CD8 + T cell response in this patient substantiallycontributed to the observed delay in recovery. Further investigation of this patient’s clinicalcourse and clinical samples is currently being undertaken.
References
1. Taubenberger, J. K., and Morens, D. M. (2008). The pathology of influenza virus infec-tions.
Annu Rev Pathol. et al. (2011). Immune responses to influenza virus infection. Virus Res. et al. (2000). Interferons: cell signalling, immune modulation, antiviralresponse and virus countermeasures.
J Gen Virol.
NatRev Immunol. et al. (1999). Natural killer cells in antiviral defense: function and regulationby innate cytokines. Annu Rev Immunol.
Annu Rev Immunol.
Nat RevImmunol. et al. (1998). Counting antigen-specific CD8 T cells: a reevaluationof bystander activation during viral infection.
Immunity.
J Virol.
Trends Immunol. + T cells: foot soldiers of the immune system.
Immunity.
J Immunol. + T cells are complementary andessential for natural resistance to a highly lethal cytopathic virus.
J Immunol.
J Exp Med. + cytotoxicT cells to virus clearance or pathologic manifestations of influenza virus infection in a Tcell receptor transgenic mouse model. J Exp Med. et al. (2010). Protective efficacy of cross-reactive CD8 + T cells recog-nising mutant viral epitopes depends on Peptide-MHC-I structural interactions and T cellactivation threshold.
PLoS Pathog. e1001039. doi:10.1371/journal.ppat.1001039.17. Yap, K., and Ada, G. (1978). Cytotoxic T cells in the lungs of mice infected with aninfluenza A virus.
Scand J Immunol. et al. (1981). Recovery from a viral respiratory infection: I. influenza pneu-monia in normal and T-deficient mice.
J Immunol. et al. (1992). Delayed clearance of sendai virus in mice lacking class I MHC-restricted CD8 + T cells.
J Immunol. et al. (1992). Transgenic mice lacking class I major histocompatibilitycomplex-restricted T cells have delayed viral clearance and increased mortality after in-fluenza virus challenge.
J Exp Med. et al. (2015). Recovery from severe H7N9 disease is associated with di-verse response mechanisms driven by CD8 + T-cells.
Nat Commun. ,7783; 10.1038/n-comms7833.22. Perelson, A. S. et al. (1996). HIV-1 dynamics in vivo: virion clearance rate, infected celllife-span, and viral generation time. Science.
Nat Rev Immunol. et al. (2003). Models of CD8 + responses: 1. what is the antigen-independentproliferation program. J Theor Biol. et al. (2004). A stochastic model of cytotoxic T cell responses.
J Theor Biol.
IntJ Immunopharmac. : 523–526.27. De Boer, R. J. (2007). Understanding the failure of CD8 + T-cell vaccination against simi-an/human immunodeficiency virus.
J Virol.
J Theor Biol. et al. (2007). Dynamics of CD8 + T cell responses during acute and chroniclymphocytic choriomeningitis virus infection.
J Immunol. et al. (2015). Mathematical modeling provides kinetic details of thehuman immune response to vaccination.
Front Cell Infect Microbiol. Theor Popul Biol.
J Theor Biol.
J Theor Biol. et al. (2007). A dynamical model of human immune response to influenzaA virus infection.
J Theor Biol.
J Virol. et al. (2009). Simulation and prediction of the adaptive immune response toinfluenza A virus infection.
J Virol. et al. (2010). Dynamics of influenza virus infection and pathology.
J Virol. et al. (2010). Quantifying the early immune response and adaptive immuneresponse kinetics in mice infected with influenza A virus.
J Virol. et al. (2013). Assessing mathematical models of influenzainfections using features of the immune response.
PLoS ONE. e0057088.doi:10.1371/journal.pone.0057088.40. Reperant, L. A. et al. (2014). The immune response and within-host emergence of pan-demic influenza virus.
Lancet. et al. (2015). Predicting pathogen-specific CD8 T cell immune responses froma modeling approach.
J Theor Biol. et al. (2016). Mathematical model reveals the role of memoryCD8 T cell populations in recall responses to influenza.
Front Immunol. et al. (2015). The time-interval between infections and viral hierarchies aredeterminants of viral interference following influenza virus infection in a ferret model. JInfect Dis. et al. (1988). Passive serum antibody causes temporary recovery from influenzavirus infection of the nose, trachea and lung of nude mice.
Immunol. et al. (2012). Modeling within-host dynamics of influenzavirus infection including immune responses.
PLoS Comput Biol. e1002588.doi:10.1371/journal.pcbi.1002588. 216. Cao, P. et al. (2015). Innate immunity and the inter-exposure interval determine the dy-namics of secondary influenza virus infection and explain observed viral hierarchies.
PLoS Comput Biol. e1004334. doi:10.1371/journal.pcbi.1004334.47. Baccam, P. et al. (2006). Kinetics of influenza A virus infection in humans.
J Virol. et al. (2012). Activation mechanisms of natural killer cells during influenzavirus infection.
PLoS ONE. e51858. doi:10.1371/journal.pone.0051858.49. Kaech, S. M., and Ahmed, R. (2001). Memory CD8 + T cell differentiation: initial antigenencounter triggers a developmental program in na¨ıve cells.
Nat Immunol. et al. (1999). Naive, effector, and memory CD8 T cells in protection againstpulmonary influenza virus infection: homing properties rather than initial frequencies arecrucial. J Immunol. et al. (2005). Frequency, specificity, and sites of expansion of CD8 + Tcells during primary pulmonary influenza virus infection.
J Immunol. et al. (2001). Na¨ıve CTLs require a single brief period of antigenicstimulation for clonal expansion and differentiation.
Nat Immunol. et al. (2013). Reducing uncertainty in within-host parameter estimates ofinfluenza infection by measuring both infectious and total viral load. PLoS ONE. e0064098. doi:10.1371/journal.pone.0064098.54. Smith, A. M. et al. (2010). An accurate two-phase approximate solution to an acute viralinfection model.
J Math Biol. et al. (2003). Fewer CTL, not enhanced NK cells, are sufficient forviral clearance from the lungs of immunocompromised mice.
Cellul Immunol.
J Immunol. et al. (2015). The administration of oseltamivir results in reduced effector andmemory CD8 + T cell responses to influenza and affects protective immunity.
FASEB J. pii: fj.14-260687.58. Asano, M. S., and Ahmed, R. (1996). CD8 T cell memory in B cell-deficient mice.
J ExpMed. et al. (2000). Response of naive and memory CD8 + T cells to anti-gen stimulation in vivo . Nat Immunol. + T-cell memory.
Immunol Rev. + T cells.
Cell Res.
Immunol Cell Biol.
Nat Rev Immunol. et al. (2015). Prior population immunity reduces the expected impact ofCTL-inducing vaccines for pandemic influenza control.
PLoS ONE. e0120138.doi:10.1371/journal.pone.0120138.65. Harty, J. T., and Badovinac, V. P. (2008). Shaping and reshaping CD8 + T-cell memory.
Nat Rev Immunol. et al. (1994). Virus-specific CD8 + T-cell memory determined by clonal burst.
Nature. et al. (2005). Age-related loss of na¨ıve T cells and dysregulation of T-cell/B-cell interactions in human lymph nodes.
Immunol. et al. (2007). Dramatic increase in naive T cell turnover is linked to loss ofnaive T cells from old primates.
Proc Natl Acad Sci USA. et al. (2010). Loss of naive T cells and repertoire constriction predict poorresponse to vaccination in old primates.
J Immunol. et al. (2010). Primary CTL response magnitude in mice is determinedby the extent of naive T cell recruitment and subsequent clonal expansion.
J Clin Invest. et al. (2012). Ecological analysis of antigen-specific CTL repertoires de-fines the relationship between naive and immune T-cell populations.
Proc Natl Acad SciUSA. et al. (2015). Sizing up the key determinants of the CD8 + T cell response.
Nat Rev Immunol. et al. (1986). Distinct short-lived and long-lived antibody-producing cell popula-tions.
Eur J Immunol. et al. (2015). The generation of antibody-secreting plasma cells.
Nat Rev Im-munol.
Eur J Immunol. et al. (2000). Diminished primary and secondary influenza virus-specificCD8 + T-cell responses in CD4-depleted Ig − / − mice. J Virol. et al. (2014). CD4 + T cell help guides formation of CD103 + lung-residentmemory CD8 + T cells during influenza viral infection.
Immunity. et al. (2015). Mortality amongst patients with influenza-associated se-vere acute respiratory illness, south africa, 2009-2013.
PLoS ONE. e0118884.doi:10.1371/journal.pone.0118884. 23 unding
PC is supported by an Australian Government National Health and Medical Research Council(NHMRC) funded Centre for Research Excellence in Infectious Diseases Modelling to InformPublic Health Policy. ZW is supported by an NHMRC Australia–China Exchange Fellowship.AWCY is supported by an Australian Postgraduate Award. JM and KK are support by NHMRCCareer Development and Senior Researcher Fellowships respectively. JMM is supported byan Australian Research Council Future Fellowship. The H7N9 clinical study was supportedby National Natural Science Foundation of China (NFSC) grants 81471556, 81470094 and81430030.
Acknowledgements
The authors would like to thank members of the Modelling and Simulation Unit in the Centrefor Epidemiology and Biostatistics and the School of Mathematics and Statistics, University ofMelbourne, for helpful advice on the study.
Competing interests
We have no competing interests. 24 ar. Description Value & Unit Ref. V initial viral load 10 [ u V ] estimated T initial number of epithelial cells in the URT 7 × cells (53) g T base growth rate of healthy cells 0 . − fixed p V viral production rate 210 [ u V ] cell − d − estimated p F IFN production rate 10 − [ u F ] cell − d − estimated p C naive CD8 + T cell proliferation rate 1 . − (32) p B naive B cell proliferation rate 0 .
52 d − estimated p S short-lived antibody production rate 12 [ u A ] cell − d − estimated p L long-lived antibody production rate 4 [ u A ] cell − d − estimated d V nonspecific viral clearance rate 5 d − (47) d I nonspecific death rate of infected cells 2 d − (32) d F IFN degradation rate 2 d − (45) d E death rate of effector CD8 + T cells 0.57 d − (59) d P death rate of plasma B cells 0.5 d − estimated d S short-lived antibody consumption rate 2 d − estimated d L long-lived antibody consumption rate 0.015 d − estimated b rate of viral consumption by binding to target cells 5 × − cell − d − estimated b ′ rate of infection of target cells by virus 3 × − [ u V ] − d − estimated f rate of conversion to virus-resistant state 0.33 [ u F ] − d − (45) r rate of recovery from virus-resistant state 2.6 d − (45) k S rate of viral neutralization by short-lived antibodies 0.8 [ u A ] − d − estimated k L rate of viral neutralization by long-lived antibodies 0.4 [ u A ] − d − estimated k N killing rate of infected cells by IFN-activated NKcells 2.5 [ u F ] − d − (45) k E killing rate of infected cells by effector CD8 + T cells 5 × − cell − d − estimated b Cn maximum rate of stimulation of naive CD8 + T cellsby virus 1 d − (9) b Bn maximum rate of stimulation of naive B cells by virus 0.03 d − estimated h C half-maximal stimulating viral load for naive CD8 + T cells 10 [ u V ] estimated h B half-maximal stimulating viral load for naive B cells 10 [ u V ] estimated t C delay for effector CD8 + T cell production 6 d (51) t B delay for plasma B cell production 4 d estimatedTable 1: Model parameter values obtained by fitting the model to experimental data. [ u V ] , [ u F ] and [ u A ] represent the units of viral load, IFN and antibodies respectively. [ u V ] and [ u A ] are EID / ml (50% egg infective dose) and pg / ml, consistent with the units of data. IFN isassumed to be a non-dimensionalized variable in the model, and therefore [ u F ] can be ignored.Some parameters are obtained from the literature and the rest are obtained by fitting the modelto experimental data in the paper of Miao et al. (38), except g T which is of minor importancewhen considering a single infection and is thus fixed to reduce uncertainty.25 upplementary Material The delay in activation of CD8 + T cells and B cells directly results in the delay of productionof effector CD8 + T cells and antibody-producing cells. By a time-shift transform t = t ′ + t C ,Eq. 7 in the main text becomes dE ( t ′ + t C ) dt ′ = b Cn ( V ( t ′ ) V ( t ′ ) + h C ) C n ( t ′ ) e ( p C t C ) − d E E ( t ′ + t C ) , (S1)where t ≥ t ′ starts from − t C . As defined in the main text, V = V ( t ′ ) = − t C < t ′ <
0. Moreover, E ( t ′ + t C ) = t C < t ′ <
0. Thus, the equation istrivial for t ′ ≥
0. Therefore, for t ′ ≥
0, we replace t ′ back to t and Eq. S1 becomes dE ( t + t C ) dt = b Cn ( V ( t ) V ( t ) + h C ) C n ( t ) e ( p C t C ) − d E E ( t + t C ) , (S2)where t ≥
0. Similarly, Eqs. 9-11 in the main text can be changed to dP ( t + t B ) dt = b Bn ( V ( t ) V ( t ) + h B ) B n ( t ) e ( p B t B ) − d P P ( t + t B ) , (S3) dA S ( t + t B ) dt = p S P ( t + t B ) − d S A S ( t + t B ) , (S4) dA L ( t + t B ) dt = p L P ( t + t B ) − d L A L ( t + t B ) , (S5)In this way, the model presented in the main text is equivalent to the following model: dVdt = p V I − d V V − k S VA S ( t ) − k L VA L ( t ) − b V T , (S6) dTdt = g T ( T + R )( − T + R + IT ) − b ′ V T + r R − f F T , (S7) dIdt = b ′ V T − d I I − k N IF − k E IE ( t ) , (S8) dFdt = p F I − d F F , (S9) dRdt = f FT − r R , (S10) dC n dt = − b Cn ( VV + h C ) C n , (S11) dE ( t + t C ) dt = b Cn ( VV + h C ) C n e ( p C t C ) − d E E ( t + t C ) , (S12) dB n dt = − b Bn ( VV + h B ) B n , (S13) dP ( t + t B ) dt = b Bn ( VV + h B ) B n e ( p B t B ) − d P P ( t + t B ) , (S14) dA S ( t + t B ) dt = p S P ( t + t B ) − d S A S ( t + t B ) , (S15) dA L ( t + t B ) dt = p L P ( t + t B ) − d L A L ( t + t B ) . (S16)26or variables whose independent variable is not explicitly specified, they are all functions of t ,i.e. V reads V ( t ) . This model avoids negative time and is solved by the following steps: • Firstly choosing a time step size D t and using it to discretise the time domain to be t = , D t , D t , D t , ..., N D t . The results in the main text are generated using D t = . D t does notimprove the solution (results not shown). • Given initial condition ( V , T , I , F , R , C n , E , B n , P , A S , A L ) = ( V , × , , , , , , , , , ) ,we solve the model iteratively. For iteration from k D t to ( k + ) D t , A S ( t ) , A L ( t ) and E ( t ) in Eqs. S6 and S8 are already known and thus treated as parameters. The systembecomes an ODE system and can be easily solved by using a built-in ODE solver ode15s in MATLAB(R2014b) with default settings. All the variables at time ( k + ) D t are thenupdated for use in the next iteration.MATLAB code is provided below. c l e a r d t = 0 . 1 ; % t i m e s t e p s i z e t i m e = 0 : d t : 1 0 0 ; % p a r a m e t e r s T0=7 e + 7 ; gT = 0 . 8 ; pV = 2 1 0 ; d e l t a V = 5 . 0 ; b e t a =5 e −
7; b e t a p =3e − d e l t a I = 2 ; kappaN = 2 . 5 ; k a p p a E =5 e − p h i = 0 . 3 3 ; r h o = 2 . 6 ; pF =1e −
5; d e l t a F = 2 ; b e t a C n = 1 ; b e t a B n = 0 . 0 3 ; k a p p a S = 0 . 8 ; k a p p a L = 0 . 4 ; pC = 1 . 2 ; pB = 0 . 5 2 ; d e l t a E = 0 . 5 7 ; d e l t a P = 0 . 5 ; pS = 1 2 ; pL = 4 ; d e l t a S = 2 ; d e l t a L = 0 . 0 1 5 ; t a u C = 6 ; t a u B = 4 ; hC=1 e + 4 ; hB=1 e + 4 ; % i n d e x i n d i c a t i n g when d e l a y e d p r o c e s s s t a r t s i n d C = r o u n d ( t a u C / d t + 1 ) ; % f o r t a u C i n d B = r o u n d ( t a u B / d t + 1 ) ; % f o r t a u B % v a r i a b l e v e c t o r s a n d i n i t i a l c o n d i t i o n s V= z e r o s ( 1 , l e n g t h ( t i m e ) ) ; V ( 1 ) =1 e + 4 ; T=V ; T ( 1 ) =7 e + 7 ; I =T ; I ( 1 ) = 0 ; R= I ; F= I ; Cn =100 ∗ o n e s ( 1 , l e n g t h ( t i m e ) ) ; Bn =100 ∗ o n e s ( 1 , l e n g t h ( t i m e ) ) ; E= z e r o s ( 1 , i n d C + l e n g t h ( t i m e ) ) ; P= I ; AS= z e r o s ( 1 , i n d B + l e n g t h ( t i m e ) ) ; AL=AS ; 27 i n i t = [V ( 1 ) , T ( 1 ) , I ( 1 ) ,R ( 1 ) , F ( 1 ) , Cn ( 1 ) , E ( 1 ) , Bn ( 1 ) , P ( 1 ) ,AS ( 1 ) ,AL( 1 ) ] ’ ; o p t i o n s = o d e s e t ( ’ R e l T o l ’ , 1 e − −
6) ; f o r i = 2 : l e n g t h ( t i m e ) [ ˜ , Y] = o d e 1 5 s ( @ODEmodel , [ 0 d t ] , i n i t , o p t i o n s , E ( i ) ,AL( i ) ,AS ( i) , tauC , tauB , p h i , r h o , d e l t a F , gT , pF , pV , b e t a , b e t a p , kappaN ,d e l t a V , d e l t a I , b e t a C n , b e t a B n , kappaE , kappaS , pL , pS , d e l t a L ,d e l t a S , d e l t a P , d e l t a E , pC , pB , kappaL , hC , hB , T0 ) ; V( i ) =Y( end , 1 ) ; T ( i ) =Y( end , 2 ) ; I ( i ) =Y( end , 3 ) ; R ( i ) =Y( end , 4 ) ; F ( i ) =Y( end , 5 ) ; Cn ( i ) =Y( end , 6 ) ; Bn ( i ) =Y( end , 8 ) ; P ( i ) =Y( end , 9 ) ; E ( i n d C + i ) =Y( end , 7 ) ; AS ( i n d B + i ) =Y( end , 1 0 ) ; AL( i n d B + i ) =Y( end , 1 1 ) ; i n i t =Y( end , : ) ’ ; % i n i t i a l c o n d i t i o n f o r n e x t i t e r a t i o n e n dIn the command of ode15s , a function “ODEmodel” is required and provided below. f u n c t i o n ynew=ODEmodel ( ˜ , y , E , AL , AS , tauC , tauB , p h i , r h o , d e l t a F , gT, pF , pV , b e t a , b e t a p , kappaN , d e l t a V , d e l t a I , b e t a C n , b e t a B n , kappaE ,kappaS , pL , pS , d e l t a L , d e l t a S , d e l t a P , d e l t a E , pC , pB , kappaL , hC , hB ,T0 ) % V : v i r a l l o a d % T : t a r g e t c e l l % I : i n f e c t e d c e l l % R : R e s i s t a n t c e l l % F : IFN % Cn : n a i v e CD8+ T c e l l s % E : e f f e c t o r CD8+ T c e l l s % Bn : n a i v e B c e l l s % P : p l a s m a B c e l l s % AS : s h o r t − l i v e d a n t i b o d i e s % AL : l o n g − l i v e d a n t i b o d i e s % y = [V , T , I , R , F , Cn , E , Bn , P , AS , AL] ynew= z e r o s ( 1 1 , 1 ) ; ynew ( 1 ) =pV ∗ y ( 3 ) − d e l t a V ∗ y ( 1 ) − k a p p a S ∗ y ( 1 ) ∗ AS − k a p p a L ∗ y ( 1 ) ∗ AL − b e t a ∗ y ( 1 ) ∗ y ( 2 ) ; ynew ( 2 ) =gT ∗ ( y ( 2 ) +y ( 4 ) ) ∗ (1 − ( y ( 2 ) +y ( 3 ) +y ( 4 ) ) / T0 ) − b e t a p ∗ y ( 1 ) ∗ y ( 2 )+ r h o ∗ y ( 4 ) − p h i ∗ y ( 2 ) ∗ y ( 5 ) ; 28 ynew ( 3 ) = b e t a p ∗ y ( 1 ) ∗ y ( 2 ) − d e l t a I ∗ y ( 3 ) − kappaN ∗ y ( 3 ) ∗ y ( 5 ) − k a p p a E ∗ y( 3 ) ∗ E ; ynew ( 4 ) = p h i ∗ y ( 2 ) ∗ y ( 5 ) − r h o ∗ y ( 4 ) ; ynew ( 5 ) =pF ∗ y ( 3 ) − d e l t a F ∗ y ( 5 ) ; ynew ( 6 ) = − b e t a C n ∗ y ( 1 ) / ( y ( 1 ) +hC ) ∗ y ( 6 ) ; ynew ( 7 ) = b e t a C n ∗ y ( 1 ) / ( y ( 1 ) +hC ) ∗ y ( 6 ) ∗ e x p ( pC ∗ t a u C ) − d e l t a E ∗ y ( 7 ) ; ynew ( 8 ) = − b e t a B n ∗ y ( 1 ) / ( y ( 1 ) +hB ) ∗ y ( 8 ) ; ynew ( 9 ) = b e t a B n ∗ y ( 1 ) / ( y ( 1 ) +hB ) ∗ y ( 8 ) ∗ e x p ( pB ∗ t a u B ) − d e l t a P ∗ y ( 9 ) ; ynew ( 1 0 ) =pS ∗ y ( 9 ) − d e l t a S ∗ y ( 1 0 ) ; ynew ( 1 1 ) =pL ∗ y ( 9 ) − d e l t a L ∗ y ( 1 1 ) ; 29 .2 Details of fitting the model to data The model contains 11 equations and 30 parameters (see Table 1 in the main text). This rep-resents a serious challenge in terms of parameter estimation, and clearly prevents a straightfor-ward application of standard statistical techniques. However, based on an extensive survey ofthe experimental literature, we have been able to identify plausible, but by no means unique,combinations of parameters that successfully explain the available data. A number of param-eters were taken directly from the literature, as per the citations in Table 1. The rest (18 pa-rameters) were estimated by calibrating the model to the published data from Miao et al. (38)who measured viral titre, CD8 + T cell counts and IgM and IgG antibodies in laboratory mice(exhibiting a full immune response) over time during primary influenza H3N2 virus infection(shown in Fig. 2 in the main text). Note that the data were presented in scatter plots in theoriginal paper (38), while we presented the data in Mean ± SD at each data collection timepoint (as shown in Fig. 2 in the main text) and fit our mean-field mathematical model to themeans.To obtain the set of parameters used for the main analysis from the experimental data in thepaper of Miao et al. (38), we took the following approach: • We first manually determined a set of the 18 parameters which produced a model solutionthat reasonably matched the experimental data shown in Fig. 2 in the main text. In detail,the main criteria include: 1) the viral load starts from about 10 EID / ml, reaches apeak of about 10 EID / ml at 1–2 days p.i., then declines rapidly from about 5 days p.i.(note that the last three data points were not considered due to the limit of detection);2) the CD8 + T cell count starts to increase rapidly at about 6 days p.i., reaches a peakof about 10 –10 at 8–10 days p.i. and returns back to zero after about 20 days p.i.; 3)the IgM level starts to increase rapidly at 4–5 days p.i., reaches a peak of about 200–300 pg/ml at about 10 days p.i. and returns back to baseline after about 20 days p.i.;4) the IgG level starts to increase rapidly at 4–5 days p.i., reaches a peak of about 800–1000 pg/ml at about 20 days p.i. and decays slowly. Given the high-dimensionality of theparameter space and limited experimental data, the procedure was essential in allowing usto identify a candidate parameter set which was not far from generating a local minimumin the following optimisation process. • The candidate parameter set was then used as an initial estimate for optimization usingMATLAB’s built-in function fmincon with default settings. The target of optimizationwas to minimize the least-squares error (LSE): L = L V + L E + L M + L G , (S17)where the four LSE components for viral load, effector CD8 + T cells, IgM and IgG weregiven respectively by L V = (cid:229) i [( V model ( t i ) − V data ( t i )) w V ( t i ) / s V ] , (S18) L E = (cid:229) i [( E model ( t i ) − E data ( t i )) w E ( t i ) / s E ] , (S19) L M = (cid:229) i [( M model ( t i ) − M data ( t i )) w M ( t i ) / s M ] , (S20) L G = (cid:229) i [( G model ( t i ) − G data ( t i )) w G ( t i ) / s G ] . (S21)30 model , E model , M model and G model indicate the model solution for variables V , E , A S and A L evaluated at time t i respectively. V data , E data , M data and G data indicate the associateddata. i is the index of the time point and the t i may differ for the four components. s V = , s E = × , s M =
300 and s G =
900 were used to scale the errors to the sameorder of magnitude. Due to the fact that data were not collected at a fixed frequency(i.e. the time interval between adjacent data points is not constant), the errors at differenttime points were assigned different weights based on the length of time intervals betweenadjacent points. For example, if viral load was measured at time t i , ( i = , , , ..., k ) , theweight function w V ( t i ) for interior points was given by w V ( t i ) = t i + − t i − ( t k − t ) , (S22)and for boundary points by w V ( t ) = t − t ( t k − t ) and w V ( t k ) = t k − t k − ( t k − t ) . (S23)This error weighting was used to weaken the domination of dense data points on modelfits. It is evident that a lot of measurements were done within the first 10 days post-infection but only a few were performed after day 20 post-infection, in particular for IgGdata (see Fig. 2 in the main text). We found that using equally weighted LSEs led to amodel fit that manifestly failed to capture those sparse data points which we believe areequally, if not more important from a more biological perspective, in determining the IgGkinetics (see Fig. S8, compared with Fig. 2 in the main text). • The parameter constraints when using fmincon were set to V ∈ [ , ] , p V ∈ [ , ] , b ∈ [ , ] , b ′ ∈ [ , ] , p F ∈ [ , ] , k S ∈ [ , ] , k L ∈ [ , ] , k E ∈ [ , ] , b Bn ∈ [ , ] , p B ∈ [ , ] , d P ∈ [ , ] , d S ∈ [ , ] , p S ∈ [ , ] , p L ∈ [ , ] , d L ∈ [ , . ] , t B ∈ [ , ] , h C ∈ [ , ] and h B ∈ [ , ] . • After obtaining a locally optimized solution for the candidate parameter set, we thenchecked the solution generated and evaluated its biological plausibility (based on thecriteria mentioned above). This step was essential as given the over-specification of themodel (in a statistical sense), it was possible for good fitting solutions to be identified byMATLAB’s optimization algorithm, which were nonetheless biologically implausible.For example, oscillatory solutions for quantities such as IgG, while providing a “good-fit” to data, were not deemed acceptable on biological grounds (see Fig. S9 for such anexample). If the optimised solution failed our (qualitative) evaluation, we returned tothe first step and redetermined a new set of parameters as new initial estimates to beoptimized using MATLAB.This entire process was repeated to arrive at the default parameter set shown in Table 1 inthe main text. To guarantee that the default parameter set was a good choice, we further ran-domly generated 10,000 sets of parameter samples near the default parameter set (within ± fmincon to search for locallyoptimized solutions. Of these, 31 generated a better LSE but all failed to meet the criteriamentioned above (results not shown).Fig. 2 in the main text shows how the model reproduces the key dynamic behavior shown inthe data. We also show in the Results section in the main text, that the model behavior is robustto perturbation of model parameters and that model predictions are reasonably consistent with31 range of other experimental data, demonstrating the plausibility, if not uniqueness (of course),of the parameter set. We emphasize that, although the default parameter set is successful in re-producing a multitude of experimental observations (e.g. full immune, knockout, re-infection)as presented in the main text, we by no means claim that this parameter set is unique. It remainsan open and challenging problem to reliably identify a biologically plausible and statisticallyidentifiable solution for what is a highly complex system, where we are severely limited byavailable experimental data. 32 .3 The model with memory
CD8 + T cells
Incorporating memory CD8 + T cells into the model in the main text, we only make twochanges. The first is adding two equations to describe the memory CD8 + T cell ( C m ) pro-liferation/differentiation, similar to Eqs. 6 and 7 in the main text dC m dt = − b Cm ( VV + h Cm ) C m , (S24) dE m dt = b Cm ( V ( t − t Cm ) V ( t − t Cm ) + h Cm ) C m ( t − t Cm ) e ( p Cm t Cm ) − d E E m . (S25)Then we change the term k E IE in Eq. 3 in the main text to k E I ( E + E m ) . Hence, similar to theapproach mentioned above that moving the delayed term from viral load to effector cells, wewrite down the model in an equivalent form, dVdt = p V I − d V V − k S VA S ( t ) − k L VA L ( t ) − b V T , (S26) dTdt = g T ( T + R )( − T + R + IT ) − b ′ V T + r R − f F T , (S27) dIdt = b ′ V T − d I I − k N IF − k E I [ E ( t ) + E m ( t )] , (S28) dFdt = p F I − d F F , (S29) dRdt = f F T − r R , (S30) dC n dt = − b Cn ( VV + h C ) C n , (S31) dE ( t + t C ) dt = b Cn ( VV + h C ) C n e ( p C t C ) − d E E ( t + t C ) , (S32) dB n dt = − b Bn ( VV + h B ) B n , (S33) dP ( t + t B ) dt = b Bn ( VV + h B ) B n e ( p B t B ) − d P P ( t + t B ) , (S34) dA S ( t + t B ) dt = p S P ( t + t B ) − d S A S ( t + t B ) , (S35) dA L ( t + t B ) dt = p L P ( t + t B ) − d L A L ( t + t B ) . (S36) dC m dt = − b Cm ( VV + h Cm ) C m , (S37) dE m ( t + t Cm ) dt = b Cm ( VV + h Cm ) C m e ( p Cm t Cm ) − d E E m ( t + t Cm ) . (S38)For variables whose independent variables are not explicitly specified, they are all functionsof t , i.e. V reads V ( t ) . Memory CD8 + T cells show a shorter delay and faster proliferation thannaive CD8 + T cells (59). The shortened delay may be caused by a shortened lag time to the firstdivision and/or a reduced delay for effector cells migrating from the lymphatic compartmentto the lung (59, 36, 42). The former reduction is about 15 hours (59) and the latter is less thanabout 12 hours (36). Thus, we choose t Cm = t C = + T cells show a33igher division rate and a lower loss rate than naive CD8 + T cells (59), based on which the netproduction rate of effector cells for memory CD8 + T cells is estimated to be about 1.5 timesof that for naive cells. Thus,we choose p Cm = . p C = . − ). In the absence of data,we assume b Cm = b C and h Cm = h C . The initial number of memory CD8 + T cells is varied asspecified in the main text or figures. We assume that the effector CD8 + T cells produced byeither naive or memory CD8 + T cells are functionally identical (i.e. then have the same decayrate d E and killing rate k E ). Note that we do not model the process of differentiation of effectorCD8 + T cells into memory cells but use a memory cell pool as an initial condition to simulateviral re-infection.MATLAB code is provided below. c l e a r d t = 0 . 1 ; % t i m e s t e p s i z e t i m e = 0 : d t : 1 0 0 ; % p a r a m e t e r s T0=7 e + 7 ; gT = 0 . 8 ; pV = 2 1 0 ; d e l t a V = 5 . 0 ; b e t a =5 e −
7; b e t a p =3e − d e l t a I = 2 ; kappaN = 2 . 5 ; k a p p a E =5 e − p h i = 0 . 3 3 ; r h o = 2 . 6 ; pF =1e −
5; d e l t a F = 2 ; b e t a C n = 1 ; b e t a B n = 0 . 0 3 ; k a p p a S = 0 . 8 ; k a p p a L = 0 . 4 ; pC = 1 . 2 ; pB = 0 . 5 2 ; d e l t a E = 0 . 5 7 ; d e l t a P = 0 . 5 ; pS = 1 2 ; pL = 4 ; d e l t a S = 2 ; d e l t a L = 0 . 0 1 5 ; t a u C = 6 ; t a u B = 4 ; hC=1 e + 4 ; hB=1 e + 4 ; betaCm = 1 ;pCm = 1 . 8 ; tauCm = 5 ; % i n d e x i n d i c a t i n g when d e l a y e d p r o c e s s s t a r t s i n d C = r o u n d ( t a u C / d t + 1 ) ; % f o r t a u C indCm= r o u n d ( tauCm / d t + 1 ) ; % f o r tauCm i n d B = r o u n d ( t a u B / d t + 1 ) ; % f o r t a u B % v a r i a b l e v e c t o r s a n d i n i t i a l c o n d i t i o n s V= z e r o s ( 1 , l e n g t h ( t i m e ) ) ; V ( 1 ) =1 e + 1 ; T=V ; T ( 1 ) =7 e + 7 ; I =T ; I ( 1 ) = 0 ; R= I ; F= I ; Cn =100 ∗ o n e s ( 1 , l e n g t h ( t i m e ) ) ; Bn =100 ∗ o n e s ( 1 , l e n g t h ( t i m e ) ) ; E= z e r o s ( 1 , i n d C + l e n g t h ( t i m e ) ) ; P= I ; AS= z e r o s ( 1 , i n d B + l e n g t h ( t i m e ) ) ; AL=AS ; Cm=5000 ∗ o n e s ( 1 , l e n g t h ( t i m e ) ) ; Em= z e r o s ( 1 , indCm+ l e n g t h ( t i m e ) ) ; i n i t = [V ( 1 ) , T ( 1 ) , I ( 1 ) ,R ( 1 ) , F ( 1 ) , Cn ( 1 ) , E ( 1 ) , Bn ( 1 ) , P ( 1 ) ,AS ( 1 ) ,AL( 1 ) , Cm( 1 ) ,Em ( 1 ) ] ’ ; o p t i o n s = o d e s e t ( ’ R e l T o l ’ , 1 e − −
6) ; f o r i = 2 : l e n g t h ( t i m e ) [ ˜ , Y] = o d e 1 5 s ( @ODEmodel with memory , [ 0 d t ] , i n i t , o p t i o n s , E34 i ) ,AL( i ) , AS ( i ) ,Em( i ) , tauC , tauB , p h i , r h o , d e l t a F , gT , pF , pV ,b e t a , b e t a p , kappaN , d e l t a V , d e l t a I , b e t a C n , b e t a B n , kappaE ,kappaS , pL , pS , d e l t a L , d e l t a S , d e l t a P , d e l t a E , pC , pB , kappaL , hC, hB , T0 , betaCm , pCm , tauCm ) ; V( i ) =Y( end , 1 ) ; T ( i ) =Y( end , 2 ) ; I ( i ) =Y( end , 3 ) ; R ( i ) =Y( end , 4 ) ; F ( i ) =Y( end , 5 ) ; Cn ( i ) =Y( end , 6 ) ; Bn ( i ) =Y( end , 8 ) ; P ( i ) =Y( end , 9 ) ; E ( i n d C + i ) =Y( end , 7 ) ; AS ( i n d B + i ) =Y( end , 1 0 ) ; AL( i n d B + i ) =Y( end , 1 1 ) ; Cm( i ) =Y( end , 1 2 ) ; Em( indCm+ i ) =Y( end , 1 3 ) ; i n i t =Y( end , : ) ’ ; % i n i t i a l c o n d i t i o n f o r n e x t i t e r a t i o n e n dThe function “ODEmodel with memory” is provided below. f u n c t i o n ynew= ODEmodel New with memory ( ˜ , y , E , AL , AS , Em , tauC ,tauB , p h i , r h o , d e l t a F , gT , pF , pV , b e t a , b e t a p , kappaN , d e l t a V , d e l t a I, b e t a C n , b e t a B n , kappaE , kappaS , pL , pS , d e l t a L , d e l t a S , d e l t a P ,d e l t a E , pC , pB , kappaL , hC , hB , T0 , betaCm , pCm , tauCm ) % V : v i r a l l o a d % T : t a r g e t c e l l % I : i n f e c t e d c e l l % R : R e s i s t a n t c e l l % F : IFN % Cn : n a i v e CD8+ T c e l l s % E : e f f e c t o r CD8+ T c e l l s % Bn : n a i v e B c e l l s % P : p l a s m a B c e l l s % AS : s h o r t − l i v e d a n t i b o d i e s % AL : l o n g − l i v e d a n t i b o d i e s % Cm: memory CD8+ T c e l l s % Em : e f f e c t o r CD8+ T c e l l s p r o d u c e d f r o m memory c e l l s % y = [V , T , I , R , F , Cn , E , Bn , P , AS , AL , Cm, Em] ynew= z e r o s ( 1 3 , 1 ) ; ynew ( 1 ) =pV ∗ y ( 3 ) − d e l t a V ∗ y ( 1 ) − k a p p a S ∗ y ( 1 ) ∗ AS − k a p p a L ∗ y ( 1 ) ∗ AL − b e t a ∗ y ( 1 ) ∗ y ( 2 ) ; ynew ( 2 ) =gT ∗ ( y ( 2 ) +y ( 4 ) ) ∗ (1 − ( y ( 2 ) +y ( 3 ) +y ( 4 ) ) / T0 ) − b e t a p ∗ y ( 1 ) ∗ y ( 2 )+ r h o ∗ y ( 4 ) − p h i ∗ y ( 2 ) ∗ y ( 5 ) ; ynew ( 3 ) = b e t a p ∗ y ( 1 ) ∗ y ( 2 ) − d e l t a I ∗ y ( 3 ) − kappaN ∗ y ( 3 ) ∗ y ( 5 ) − k a p p a E ∗ y( 3 ) ∗ ( E+Em) ; ynew ( 4 ) = p h i ∗ y ( 2 ) ∗ y ( 5 ) − r h o ∗ y ( 4 ) ; ynew ( 5 ) =pF ∗ y ( 3 ) − d e l t a F ∗ y ( 5 ) ; ynew ( 6 ) = − b e t a C n ∗ y ( 1 ) . / ( y ( 1 ) +hC ) ∗ y ( 6 ) ; ynew ( 7 ) = b e t a C n ∗ y ( 1 ) . / ( y ( 1 ) +hC ) ∗ y ( 6 ) ∗ e x p ( pC ∗ t a u C ) − d e l t a E ∗ y ( 7 ) ; ynew ( 8 ) = − b e t a B n ∗ y ( 1 ) . / ( y ( 1 ) +hB ) ∗ y ( 8 ) ; ynew ( 9 ) = b e t a B n ∗ y ( 1 ) . / ( y ( 1 ) +hB ) ∗ y ( 8 ) ∗ e x p ( pB ∗ t a u B ) − d e l t a P ∗ y ( 9 ) ;35 ynew ( 1 0 ) =pS ∗ y ( 9 ) − d e l t a S ∗ y ( 1 0 ) ; ynew ( 1 1 ) =pL ∗ y ( 9 ) − d e l t a L ∗ y ( 1 1 ) ; ynew ( 1 2 ) = − betaCm ∗ y ( 1 ) . / ( y ( 1 ) +hC ) ∗ y ( 1 2 ) ; ynew ( 1 3 ) = betaCm ∗ y ( 1 ) . / ( y ( 1 ) +hC ) ∗ y ( 1 2 ) ∗ e x p ( pCm ∗ tauCm ) − d e l t a E ∗ y( 1 3 ) ; 36 ime (days) f r a c t i on o f un i n f e c t ed ep i t he li a l c e ll s FIGURE S1: Model solution with parameters given in Table 1 in the main text shows that the loss ofuninfected epithelial cells is maintained within 10–20% of total cell pool. This number of uninfectedepithelial cells equals the sum of target cells T and resistant cells R shown in Fig. 3 in the main text. ime (days)0 2 4 6 8 10 l o g ( V ) ( E I D / m l ) viral loadtarget cell number time (days) l o g ( V ) ( E I D / m l ) viral loadtarget cell number (a) (b) FIGURE S2: Inclusion of IFN prevents target cell depletion. In the absence of adaptive immunity(letting C n = B n = p F = − ) simulates a sustainedelevation of viral load and large steady state target cell number (dashed curve in panel (a)), both ofwhich are consistent with previous experimental data (shown in Fig. 5 in the main text). However, in theabsence of an innate response (letting p F = (target cell number). ‡ ·
13 14 1505100 1 (cid:181) ¶ • ‚
13 14 15 l o g ( „ ) ( E I D / m l ) ” » … ‰
13 14 150510 time (days) (cid:190) ¿ (cid:192) `
13 14 15 ´ o ˆ o ˜ e v e n t s ¯ no immunity ˘ ull immuneIF ˙ only FIGURE S3: The three-phase model behavior is robust to the change of parameters. By allowing allthe parameters to vary by ±
20% from their default values (i.e. those shown in Table 1), 1000 samplesof parameter sets were selected using Latin hypercube sampling and 1000 corresponding viral loadsolutions are shown in the upper panel. In Fig. 4 in the main text, we illustrated the phase separationusing area shading. Here, the rough times of phase separation are indicated by the histograms in thelower panel; the left one indicates the time separating the first and the second phases (i.e. the time whenthe corresponding solutions of “full immune” and “no immunity” differ by 0.2 (in log-scale) for the firsttime) and right one indicates the time separating the second and the third phases (i.e. the time when thecorresponding solutions of “full immune” and “IFN only” differ by 0.2 (in log-scale) for the first time).The solutions of “no immunity” are obtained by letting p F = b Cn = b Bn = b Cn = b Bn = ime (days)0 5 10 15 20 l o g ( V ) ( E I D / m l ) C n ¨ (cid:201) n ˚ ¸ n (cid:204) ˝ n ˛ n = 0 time (days)0 5 10 15 20 l o g ( V ) ( E I D / m l ) p S = 12p S = 8p S = 6p S = 2p S = 0 time (days)0 5 10 15 20 l o g ( V ) ( E I D / m l ) p L = 4p L = 1.6p L = 0.8p L = 0.4p L = 0 (a) (b)(c) FIGURE S4: Dependence of model behaviour on the production rate of CD8 + T cells (panel (a)),short-lived antibodies (panel (b)) and long-lived antibodies (panel (c)). Different levels of productionfor CD8 + T cells, short-lived and long-lived antibodies are modelled respectively by varying the naiveCD8 + T cell number ( C n ), short-lived antibody production rate ( p S ) and long-lived antibody productionrate ( p L ). For each panel, all other model parameters were kept fixed at the values given in Table 1 inthe main text. C (days)4 5 6 — a (cid:209) e (cid:210) a g e (cid:211) D (cid:212) (cid:213)(cid:214) e ll s ( (cid:216) (cid:217) e ll s ) model sim (cid:218) latio (cid:219) e (cid:220) po (cid:221) e (cid:222) tial (cid:223) it (cid:224)Æ (days)4 5 6 (cid:226) ª e (cid:228) o (cid:229) e (cid:230) y t i m e ( d a y s ) (cid:231) Ł model sim Ø latio Œ pie º e (cid:236) ise li (cid:237) ea (cid:238) (cid:239) it a (cid:240) e æ age (cid:242) D8 (cid:243)(cid:244) ı ells ( (cid:246) (cid:247) ells)0 2 4 6 8 10 ł e ø o œ e ß y t i m e ( d a y s ) model sim (cid:252) latio (cid:253) ( (cid:254) a (cid:255) yi n g (cid:1)C ) e x po (cid:0) e (cid:2) tial f it (a) (b)(c) FIGURE S5: Dependence of the average effector CD8 + T cell number (over days 6–20) on the delayinduced by naive CD8 + T cell activation and differentiation ( t C ). Recovery time is defined to be thetime when viral load falls to 1 EID / ml. Panel (a) shows that the average effector CD8 + T cell numberis exponentially related to the delay t C . Panel (b) shows that the recovery time is related to the delay t C in an approximately piecewise linear manner. Panel (c) shows varying delay t C also preserves theapproximately exponential relationship between the average CD8 + T cell numbers and recovery time,consistent with the results of varying initial naive CD8 + T cell numbers shown in Fig. 6 in the maintext. p V = 50 p V = 100 p V = 150 r e c o v e r y t i m e ( da ys ) p V = 200 p V = 250 p V = 300 average CD8 + T cells ( × cells)0 2101520 p V = 350 p V = 400 FIGURE S6: The exponential relationship between recovery time and the average effector CD8 + T cellnumber (over days 6–20) is robust to the changes in the viral production rate p V . Recovery time isdefined to be the time when the viral load falls to 1 EID / ml. For each p V , the relationship betweenaverage CD8 + T cells and recovery time was obtained by varying initial naive CD8 + T cell number.Except for large p V (e.g. p V = , p V is large (the cases of p V = , + T cells. For example adecrease in recovery time is observed for a lower level of effector CD8 + T cells (in the range of lessthan 5000 cells) for p V = F = 0 0 2810121416 p F = 1e-9 0 2810121416 p F = 1e-80 2 r e c o v e r y t i m e ( da ys ) F = 1e-7 0 2810121416 p F = 1e-6 0 2810121416 p F = 1e-50 2810121416 p F = 1e-4 average CD8 + T cells ( × cells) F = 1e-3 FIGURE S7: The exponential relationship between recovery time and the average effector CD8 + T cellnumber (over days 6–20) is robust to the change of IFN production rate p F . Recovery time is defined tobe the time when viral load falls to 1 EID / ml. The thin black curves through the simulated data areexponential fits. One may observe that the recovery time increases and then decreases as p F increases.The reason is not entirely clear. Our explanation is that the increase for small pF is likely becauseincreasing p F reduces the extent of target cell depletion and thus makes infection longer. When pF issufficiently large, a further increase in p F has little effect on target cell pool but induces stronger innateimmune response, which play a role in shortening the recovery time.
10 20 l og ( V ) ( E I D / m l ) -202468 0 10 20 30 CD + T c e ll s × I g M ( pg / m l ) time (days) I g G ( pg / m l )05001000
10 20 l og ( V ) ( E I D / m l ) -202468 0 10 20 30 CD + T c e ll s × I g M ( pg / m l ) time (days) I g G ( pg / m l )05001000 FIGURE S8: An example of a best-fit solution failing to capture the sparse part of the IgG data verywell. Black dots are experimental data from (38) and black curves are the fits. The fit was generatedby using equally weighted LSE function (scale factors remain unchanged) and the parameters in Table1 in the main text as an initial guess in MATLAB’s fmincon function. Compared to using weighted LSE(which generates the fit shown by green dashed curves), using equally weighted LSE does not alter thefits to the viral load and CD8 + T cell number but significant underestimates antibody levels. Note thatdue to the limit of detection for the viral load (occurring after 10 days post-infection as seen in viralload data), the last three data points in the upper-left panel were not taken into consideration for modelfitting.
10 20 l og ( V ) ( E I D / m l ) -202468 0 10 20 30 CD + T c e ll s × I g M ( pg / m l ) time (days) I g G ( pg / m l )05001000