In January 2020, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) was identified as the infectious agent causing an outbreak of viral pneumonia in Wuhan, China. It was soon established that droplet-based human to human transmission was occurring, and on March 11, 2020 , the World Health Organization characterized coronavirus disease 2019 (COVID-19) as a pandemic. As of this article's submission date, COVID-19 has infected more than 44.9 million people, causing more than one million deaths. A pandemic-scale outbreak creates tremendous socioeconomic burden due to thwarted productivity, a spike in healthcare expenses, and irreparable loss of human lives 1, 2 . Furthermore, implementation of social and physical isolation measures has caused many countries to declare states of emergency and lockdowns with border closures.
SARS-CoV-2 is the seventh identified human coronavirus and the third novel one to emerge in the last 20 years. It is a single-stranded positive sense RNA genome of about 30,000 nucleotides that encodes ~27 proteins and four structural proteins. A surface-expressed spike protein mediates receptor binding and membrane fusion with host cells, and the virus interacts with the angiotensin converting enzyme 2 (ACE2) receptor to gain entry into cells 3 . ACE2 mRNA is present in almost all human organs, but the receptor is particularly highly expressed on the surface of lung alveolar epithelial cells and enterocytes of the small intestine, thereby allowing a preferential accumulation of the virus in these organs 4 . The incubation period of SARS-CoV-2 ranges from about 3-17 days, and COVID-19 diagnosis cannot be made based on symptoms alone as, most are nonspecific and may be confused for more common ailments. The more serious sequelae of infection includes acute respiratory distress syndrome (ARDS) and sepsis caused by the cytokine storm from the immune response to infection, which is believed to be the leading cause of mortality in COVID-19 patients 5 . Screening for COVID-19 is done via nucleic acid testing by RT-PCR (specimens from both upper and lower respiratory tracts) and pulmonary CT scans. The viral load in naso-or oro-pharyngeal swabs is the key clinical biomarker of COVID-19 and also the key clinical endpoint of pharmacological intervention.
Although several antiviral and immunomodulatory drugs are being used for symptomatic treatment and viral load reduction, there are still no proven therapeutics for COVID-19 to date. To explore novel and effective therapeutic targets, we require a better understanding of the pathogenesis of COVID-19, particularly of virus-host interactions 6 . This will also enable more efficient disease management strategies, such as deriving prognostic information from viral load kinetics, and quantification of the effects of the immune system in controlling the disease. With limited studies on the in vivo dynamics of SARS-CoV-2, a mathematical modeling approach can be an excellent, complementary tool for investigating viral-host interactions and simulating pathogenesis in order to better understand disease progression and evaluate treatment strategies.
Indeed, the application of mathematical modeling and quantitative methods has been instrumental in our understanding of viral-host interactions of various viruses, including influenza, HIV, HBV, and HCV 7 . These kinetic models have been developed for various spatial scales, including molecular, cellular, multicellular, organ, and organism. By analyzing viral load kinetics, these models have deepened our understanding of the fundamentals of virus-host interaction dynamics, innate and acquired immunity, mechanisms of action of drugs, and drug resistance [8] [9] [10] [11] [12] . 21 . Although also a target cell-limited model, by only including upper and lower respiratory tract compartments, this model omits key biological mechanisms involved in the complete immune response, and is thus unable to provide deeper insights into the system-wide dynamics and interplay of disease response.
In order to improve upon the existing models, we have developed a multiscale semi-mechanistic model of viral dynamics, which, in addition to capturing virus-host interactions locally, is also capable of simulating the whole-body dynamics of SARS-CoV-2 infection, and is thereby capable of providing insights into disease pathophysiology and the typical and atypical presentations of COVID-19. Importantly, using our modeling platform, we can identify treatment strategies for effective viral load suppression under various clinically relevant scenarios. We note that while the modeling platform is developed for SARS-CoV-2, we also expect it to be applicable to other viruses that have shared similarities in mechanisms of infection and physical dimensions.
We have developed a semi-mechanistic mathematical model to simulate the whole-body biodistribution kinetics of SARS-CoV-2 following infection through the nasal route (Figure 1 ;
Methods: Model development). The model was formulated as a system of ordinary differential equations (Equations 1-40) that describe cellular-scale viral dynamics, whole-body transport and excretion of viruses, and innate and adaptive immune response to predict the viral load kinetics of SARS-CoV-2 in the respiratory tract, plasma, and other organs of the body. SARS-CoV-2 exhibits ACE2 tropism 22 , therefore the organs included in the model were chosen based on the presence of ACE2 receptor expressing cells in their tissues [23] [24] [25] . Specifically, the key processes described by the model include infection of ACE-2 expressing susceptible cells by SARS-CoV-2 (also referred to as target cells), production of new virions by infected cells, death of infected cells due to cytopathic effect, transport of virions from the site of infection to other organs of the body, hepatobiliary excretion of the virions, and key processes in the innate and adaptive immune response against the virus and infected cells to clear the infection. Note that in the absence of a thorough understanding of the mechanistic underpinnings of viral shedding in the feces 26 , and a growing evidence of liver damage in COVID-19 patients 27, 28 , we assumed bile production rate as the rate limiting step in the hepatobiliary excretion of the virus into the feces.
While some of the parameters of the model were known a priori (Table 1) , the remaining parameters were estimated through nonlinear regression using published in vivo 29 and clinical 30 data. Specifically, from published experimental data for hamsters 29 , we first calibrated a reduced version of the model (referred to as Reduced model; Equations 1-23) that comprises all compartments and interferon (IFN)-mediated innate immunity, but lacks adaptive immunity (bottom half of Figure 1 ; also see workflow in Figure 2 ). The parameters of the reduced model characterize cellular-scale viral dynamics, IFN-mediated immunity, inter-compartment viral transport, and hepatobiliary excretion of the virus from the mononuclear phagocytic system (MPS). The estimated parameters were then used in the complete version of the model (referred to as Full model; Equations 1-40), which also includes adaptive immunity, to calibrate the remaining parameters using nonlinear regression with clinical data 30 . The models were solved numerically in MATLAB as an initial value problem, using the built-in stiff ODE solver ode15s.
Calibrating parameters of the reduced model As shown in Figure 3 , the numerical solution of the reduced model for whole-body viral kinetics, IFN kinetics, and target cell population kinetics in hamsters satisfies the initial conditions, and is in good agreement with the available in vivo data 29 for viral and IFN kinetics (Pearson correlation coefficient R between experimental data and model fits is > 0.97, p < 0.0001, Figure S1a) . The corresponding parameter estimates are given in Table 2 . Based on findings in the in vivo study by Chan et al. 29 that the adaptive immune response in test animals was not triggered during the first seven days post-infection, it is reasonable to use the reduced form of the model to estimate the unknown parameters, rather than using the full model at this point.
The model solution (Figure 3) shows the kinetics of ACE2-expressing target cells (solid orange lines) and their infected counterparts (dashed orange lines) in every compartment. These infected cells can produce new virions that will in turn infect other healthy target cells. Because we are . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
The copyright holder for this preprint this version posted November 3, 2020.
using a target-cell limited modeling assumption 9, 18 , the healthy target cells that become infected by the virus are not replaced by new healthy cells, and as seen in Figure 3 , the target cells were observed to deplete within 48 h post infection. The viral load kinetics (blue curve) is primarily governed by the interplay of new virion production, distribution of the virions between compartments, viral elimination by alveolar and MPS macrophages, hepatobiliary excretion of viruses from the body, cytopathic death of infected cells, and suppression of viral production due to IFN produced by infected cells 9 , which is shown in Figure 3i . As the infected cell population tapers, the IFN concentration will also decrease to the pre-infection baseline value. In our model, infected cells of the respiratory tract are the source of IFN following infection, the lack of which has been found to be the underlying cause of life-threatening COVID-19 due to uncontrolled viral replication in the absence of IFN regulation [31] [32] [33] .
Of note, the plasma compartment (Figure 3h ) of the model does not contain any target cell population and thus its viral load kinetics is only governed by the influx and outflux of viruses from various compartments. However, in the full model, the neutralization of viruses by antibodies will also be considered in the plasma compartment, as discussed in the next section. Plasma flow is the key mechanism of viral transport and systemic spread of infection in the body 34 , but due to lack of established mechanistic underpinnings of these processes, we instead use phenomenological rate constants to characterize viral transport. Based on the estimated characteristic times (1-24 h) of the vascular transport processes (shown in Figure 1 and presented as rates in Table 2 ), it can be inferred that viral transport is permeability-limited and not perfusionlimited, i.e., capillary permeability and vascular surface area govern the rate of extravasation of virions from blood vessels into tissue interstitium to reach the target cells, and thus viral transport is not exclusively governed by the plasma flow rates into the organs. This is consistent with the in vivo behavior of nanomaterials of comparable size [35] [36] [37] [38] [39] [40] , and is in contrast to the perfusion ratelimited kinetics of smaller lipophilic molecules. The variability in characteristic times of vascular transport can be explained by differences in the permeability of capillary endothelium due to differences in pore sizes of endothelial fenestrae 41 . For instance, the blood brain barrier seems to resist transport of virus to the brain, thereby leading to an estimated characteristic time of influx of 1 day, which is ~twenty-times longer than the estimated characteristic time of influx to the MPS (1.25 h) that contains large sinusoidal pores in its microvasculature. Of note, the non-vascular transport processes have relatively longer characteristic times that can be attributed to resistance to transport offered by mucus or degradation caused by pH, among other factors.
Once the parameters discussed in the previous section were estimated, they were then used in the full model to calibrate the remaining parameters (see Table 3 ) relevant to the innate and adaptive immune system using published clinical data (n = 4 untreated patients) 30 . Due to the uncertainty associated with the duration between day of infection and onset of symptoms (referred to as incubation period), a shifting parameter was included in the calibration routine. Numerically, the time points corresponding to the data were shifted units of time. As shown in Figure 4 , the model correctly represents the initial conditions of the variables and predicts an incubation period of ~6 days (indicated by red arrow in Figure 4a ), which is comparable to published literature 42 . Also, assuming nasal route as the route of infection, the numerical solution for URT at time = 0 suggests exposure to a viral load of ~10 7 copies/mL. The clinical data shows the viral load kinetics in the upper (Figure 4a ) and lower respiratory tract (Figure 4b) , the IFN kinetics (Figure 4i) , the . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
The copyright holder for this preprint this version posted November 3, 2020. ; effector CD8 + (CD8 * ) and activated CD4 + (CD4 * ) cell population kinetics (Figure 4k) , and the total neutralizing antibody kinetics (Figure 4l) . As shown in Figure 4 , the full model solution fits the data well (Pearson correlation coefficient R > 0.98, p < 0.0001, Figure S1b) , and was able to predict the kinetics of viral load in the remaining compartments by using the viral dynamics and transport parameters estimated from the in vivo data (through calibration of the reduced model).
The model predicts that the viral load in extrapulmonary organs and plasma persists for ~17-20 days post onset of symptoms, consistent with published studies 43 , and is thus comparable to the duration of viral detection in URT and LRT. Therefore, it can be also inferred that nasopharyngeal swabs can safely provide an indication of the infection status of the patients.
The model also shows the kinetics of naïve lymphocytes (Figure 4j ) and antibody producing plasma cells (Figure 4k ) in the lymphatic compartment, which is represented as a common compartment for the entire body. Importantly, in close agreement with published literature 44, 45 , the model predicts that the systemic concentration of antibodies persists above the detectable limit for >100 days post onset of symptoms, following which it may no longer be detectable ( Figure   4l ). This finding suggests the lack of indefinite antibody protection against reinfection 46 , and highlights the need for vaccine boosters to achieve long-lasting immunity 47 . We note that the data used in the above calibrations was obtained under conditions where neither the animals nor the patients were given any pharmacological treatment. Hence, the data are appropriate to calibrate the effects of the immune components and other physiological processes in distributing and eliminating the viral load.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
While URT and LRT are the preferred sites to detect the presence of SARS-CoV-2, it is important to note, and as is evident from the model predictions, that the viral load in non-pulmonary organs can attain comparable levels, and can thus explain the non-respiratory symptoms observed in some COVID-19 patients 28, 48, 49 . Following transport of the virus from respiratory tract to blood or via gastrointestinal tract to blood, organs that have a significant population of ACE2 expressing cells may become infected by the virus, leading to the extrapulmonary manifestations of COVID-19 that may include symptoms such as diarrhea and impaired renal-, hepatic-, cardiovascular-, or neurological-functions.
To investigate the effects of the cellular and humoral arms of innate and adaptive immunity on viral load kinetics in the body, we individually switched off these components and simulated the whole-body viral kinetics for up to the time when viral load fell below the detectable limit 18 of 10 2 copies ml -1 . This numerical experiment is meant to mimic the effect of compromised immunity due to an underlying condition in a virtual patient undergoing no antiviral treatment. As seen in the viral kinetics in Figure 5 , when one or all of the immune components were shut down, the viral concentration in all the compartments was higher than the baseline (dashed dark red line).
Further, it can be inferred that IFN is the primary mechanism of controlling viral load in the URT ( Figure 5a ) and GI (Figure 5c ), while macrophages (alveolar macrophages, Kupffer cells, and splenic macrophages) play a predominant role in limiting infection in the LRT, MPS, plasma, and other organs, followed by IFNs (Figures 5b, 5d-5h) . Findings in the literature support the above observations as follows. Lack of IFNs can lead to excessive viral production 33 and cause lifethreating COVID-19 in patients deficient in functional IFNs due to, for example, the occurrence . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint of loss-of-function mutations in genes governing IFN-mediated immunity 32 or auto-antibodies against IFNs 31 . Further, FABP4+ alveolar macrophages were observed to be largely absent in patients with severe COVID-19, but were a predominant macrophage in patients with mild disease 50 , indicating the major role of FABP4+ alveolar macrophages in controlling infection as also shown previously for patients with chronic obstructive pulmonary disease (COPD) 51 . We assume a constant supply of macrophages in the lungs and MPS in our model, but a deleterious effect of infection on these immune cells cannot be ruled out, and further experimental evidence is necessary to model the cell population kinetics appropriately 52 . Further, the model reveals that the effect of adaptive immunity (antibodies and CD8 * cells) is not significant in controlling infection, but it does not necessarily rule out the therapeutic potential of exogenously administered antibodies or novel cell-based therapies (e.g. T cell therapy).
All the above scenarios abstractly represent real-world underlying conditions (e.g. cancer, diabetes, autoimmune diseases) in patients that lead to varying degrees of immunosuppression, thus highlighting the importance of an individual's immune status in regulating viral kinetics. In the absence of immune responses, the only plausible mechanism that brings the viral load down is hepatobiliary excretion of the virus through the MPS, but our results show it takes several weeks before the viral load falls below the detectable value of 10 2 copies ml -1 . Note that upon shutting down IFN-and macrophage-mediated immunity or total immunity, the viral load grew beyond clinically observed values in the literature (~10 12 copies/ml), therefore we set an upper bound at 10 12 copies/ml to keep the results clinically meaningful.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. Treatment over a combination of the two parameters, i.e., IFN and viral production rates, can be more efficient than monotherapy in suppressing viral load, as also suggested by clinical studies 54, 55 . While cytopathic death rate and viral production rate play major roles in governing the clearance time of viruses ( clear ), IFN production rate is relatively less significant. Further, transport processes play a significant role in governing both total viral load and time to clear the load.
Specifically, the transport of virus from URT to LRT ( 1 ) affects viral load in URT, and blood to MPS via the hepatic artery ( ) and MPS to blood via the hepatic vein ( ) govern the clearance time from plasma, which should also affect the viral clearance of other organs connected to plasma.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. are important parameters in governing the kinetics of viral load in LRT and plasma, respectively.
Finally, the total viral load in LRT is also affected by pulmonary absorption ( 1 ), i.e. by the rate of transport of viruses from LRT to plasma.
Up to this point, we have demonstrated that the model was able to reproduce viral kinetics that were consistent with experimentally measured values and produce observations consistent with published literature. We next sought to examine whether this tool can make predictions that might allow clinicians to design an effective therapy for patients, helping to optimize and personalize their treatment regimen. As a numerical experiment, we tested three treatment scenarios in controlling infection: a hypothetical antiviral agent, interferon therapy, and antiviral agentinterferon combination therapy. Further, the effects of the timing of therapy initiation ( ) were also tested, i.e., starting therapy on the day of onset of symptoms ( = onset ), and starting therapy 5 days post onset of symptoms ( = onset + 5). To quantify treatment effectiveness, we compared the viral load kinetics in three compartments, namely, URT, LRT, and plasma in simulations spanning a period of four weeks.
We simulated therapy with a hypothetical antiviral agent that has the same mechanism of action as Remdesivir 56 , i.e. interference with RNA-dependent RNA polymerase. For this, 200 mg loading dose (≡initial plasma concentration 0 = 5 M, assuming a molecular weight of 800 g/mol and a volume of distribution of 50 L), followed by 100 mg daily maintenance doses for 9 additional . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint days, via intravenous injection were simulated to mimic the pharmacological intervention.
Treatment was started at one of the two time points previously mentioned. As a simplification, assuming one-compartment pharmacokinetics for the hypothetical drug with elimination rate constant Cl = 1 −1 , we use the plasma concentration kinetics as a surrogate for tissue concentration kinetics of the drug. The plasma concentration ( ) is thus given by,
Hence, the injection times are = , + 1, … , + 9, and we define this set of times as . The antiviral agent acts by inhibiting the production rate of virus from infected cells in the body, 
We simulated treatment with interferon beta-1a, with and without the hypothetical antiviral agent (discussed above). For this, four subcutaneous injections (Dose = 44 µg each) of interferon beta-. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ;  1a were administered every other day starting at one of the two time points previously mentioned.
The plasma concentration kinetics of interferon beta-1a I ( ) is obtained by solving the following equation:
Hence, the injection times are = , + 2, … , + 6, and we define this set of times as . The . Figure 7d -f, interferon beta-1a therapy significantly reduces the viral load, and the effect is comparable to the hypothetical antiviral agent. Also, an early initiation of therapy leads to a lesser total viral load and a faster reduction in viral load.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint Combination therapy Finally, we tested the effect of combination therapy with the hypothetical antiviral agent and interferon beta-1a on viral load kinetics (plasma concentration kinetics in the inset of Figure 7g ). Figure 7g -i, the combination therapy leads to a faster reduction in viral load compared to the effect of antiviral agent and interferon beta-1a alone, such that the viral load seems to fall below the detection threshold 2-3 days earlier.
In summary, we have developed a semi-mechanistic mathematical model that predicts the whole- Through global sensitivity analysis, we identified the key mechanisms that control infection and can be used as potential therapeutic targets for pharmacological intervention. Finally, we tested the potential of such therapeutic targets by simulating clinically relevant treatment options and identified the importance of the timing of treatment initiation, and the effects of various therapies to suppress infection effectively and immediately. As a limitation of the model, it is important to note that due to lack of clinical data for the viral load kinetics in extrapulmonary compartments, we had to rely on the transport parameters estimated from in vivo data, and the model predictions for the extrapulmonary compartments could not be clinically validated. As clinical knowledge of the disease and its mechanisms improves, we will continue to fine-tune the model and integrate it . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint with a physiologically-based pharmacokinetic model to build an in-silico platform that can be used to simulate disease progression and a more complete pharmacokinetics of various test drugs and novel formulations.
In each organ compartment, we consider the kinetics of three populations. The first is the . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint The rate of change of total infected cell population within an organ is a function of three mechanisms. The first is an increase due to freshly infected healthy cells by viral particles, the second is death due to cytopathic effect of viral infection, proportional to the number of infected cells and characterized by the rate constant 1 , and the third is due to the interaction between effector CD8 + cells (concentration denoted by CD8 * ( )) and the infected cells. The cytotoxic effect is measured in terms of the death rate constant CD8 . It is important to note that in the reduced form of the model, infected cell death due to effector CD8 * cells is not included. In our model description, all activated immune cells are indicated with an asterisk (e.g., CD8 * ), while inactive or Naïve populations are indicated by standard naming conventions (e.g., CD8 + ).
The concentration of viral particles in each organ compartment is influenced by different factors, some of which are organ-specific. However, in all organ compartments, the rate of change of viral particles is proportional to the number of infected cells with virus production constant , but is inversely proportional to the concentration of IFN, denoted IFN( ). IFN is part of the innate immune response that acts by suppressing the viral production rate of infected cells, and is controlled by the effectiveness constant . In all compartments, viral particles may be neutralized due to the aggregation of antibodies on their surface proteins. The antibody concentration is represented by Ab( ), and the destruction rate of viral particles due to antibodies is denoted by Ab . Of note, neutralization of virus by antibodies is not incorporated in the reduced form of the model. Lastly, the viral load can also be reduced by tissue resident macrophages in the lungs (alveolar macrophages) and MPS (Kupffer cells and splenic macrophages) that engulf the viral particles and destroy them. The rate of removal of viral particles by macrophages is given by ϕ .
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint We introduce the system of ordinary differential equations that governs the concentration kinetics as follows. Mechanisms particular to each compartment are discussed after each set of equations.
Equations for the URT: which is quantified through the pulmonary absorption coefficient 1 .
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. Additionally, as discussed above, viral particles from the URT are transported to the GI tract at a rate 2 .
Equations for the MPS: is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. Given that the size of the SARS-CoV-2 virus is ~100 nm, we assume that renal excretion is not feasible 58 . Therefore, the viral kinetics in kidneys is dependent on the standard mechanisms and is affected by influx via the renal arteries ( ) and outflux via the renal veins ( ).
Equations for the Brain: While the blood brain barrier may be deterrent to the establishment of infection in the brain, we have included the brain compartment due to the lack of evidence for the former 59 is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint The equation for the plasma compartment incorporates all the outgoing fluxes of viral particles through arteries ( , , CO , CA ) and the incoming fluxes via veins ( , , CO , ). Also, the viral load from LRT in equation (6) gets added to plasma at rate 1 . Lastly, similar to all the previous compartments, viral particles can be neutralized by antibodies at a rate Ab .
The model up to this point is the reduced model, and the only equation that captures the immune system is equation (7). This equation describes the change of concentration of IFN, which is part of the innate immune system. However, the adaptive immune system should activate to properly mount a full immune response. The equations that follow provide the remaining elements to initiate and maintain a humoral and cell-mediated adaptive immune response to make up the full model. if the difference is large, the denominator is equally large, and the ratio is close to one. This allows a near constant modulation rate that only acts whenever the population of macrophages is far from equilibrium. This reequilibration rate is denoted by APC . The APC population is also impacted by the interaction between APCs and invading viral particles. This causes the population of APCs to . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ;
Furthermore, each of these lymphocytes can be promoted to its activated state, which we denote with an asterisk, * . This conversion takes place when activated APCs interact with naïve lymphocytes, presenting them fragments of the ingested viral particles. The transformation rate from the naïve state to the activated one is represented by the constant , where indicates the cell type. For the case of CD4 + and CD8 + , we use the common rate T . The concentration of activated lymphocytes can fluctuate due to two mechanisms. One is the conversion of naïve cells, which was already discussed. The second is due to death caused by a variety of factors like exhaustion or apoptosis, and is expressed through the death constant . B cells, however, follow a different mechanism, which we discuss later. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint The final equation in the model describes the rate of change of antibody concentration, Ab( ). The sole contributors to the production of antibodies are short-lived and long-lived plasma cells. The corresponding production rates are given by P S and P L , respectively. We consider two mechanisms by which antibodies can be consumed or eliminated. One is the loss of antibodies due to their interaction with viral particles in a given compartment with rate constant Ab . The second is the elimination of antibodies by other factors independent of the viral load, e.g., degradation or clearance from tissues, which occurs at rate constant Cl Ab . A summary of the parameters used in the model is given in Tables 1-3.
Once the full set of parameters for the model has been defined and fixed for the conditions corresponding to Figure 4 , we proceed to identify the relative effect of each parameter on the viral kinetics. This is done by systematically perturbing one or multiple parameters and comparing the outcomes between the perturbed state and the original state. To quantify the differences between the two states we introduce the following 6 criteria. The first three refer to the area under the curve The next three denote the total time required for the viral load in the URT, LRT, and Plasma compartment to reduce to < 10 2 viral copies per mL, i.e., 2 in log base 10 units. The last comparison criterion is the antibody titer at day 30 in the body, i.e., Ab( = 30). In the next two . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint subsections we compare each of the 6 criteria using global sensitivity analysis (GSA), applied to the 31 parameters.
Local sensitivity analysis only allows one parameter to change at a time, and thus may not be suitable for analyzing complex biological systems. We performed global sensitivity analysis (GSA 36, 54, 55 ; i.e., multiple parameters can be varied simultaneously) to obtain a further understanding of the relationship between parameters and their combined effects on the model outputs of interest. Specifically, in our analysis, all parameters were perturbed at the same time, and we assume a uniform distribution for each parameter. The range of values of each sample was PR ± 99%(PR), where PR is the reference value of a given parameter. To ensure a complete sampling of the full parameter space we used a Latin hypercube sampling scheme. Subsequently, each sample was used to recompute the kinetics, from which we extract the 6 criteria mentioned above. In total, 10 batches of 5000 samples each were computed. Then, the data were subjected to a multivariate linear regression analysis to produce an impact factor, which we label as the sensitivity index (SI). Lastly, the sensitivity indices were ranked by means of a one-way ANOVA and a Tukey's test. The higher the SI, the more significance a parameter holds. A detailed description of the GSA workflow can be found in 38, 60, 61 . All analyses were performed in MATLAB R2018a.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint Figure 2 . Modeling workflow. In the first step, the reduced model uses the hamster data to fit parameters corresponding to viral dynamics, transport coefficients, and innate immunity.
Subsequently, the more complex full model uses the previously estimated parameters and clinical data to estimate the remaining parameters pertaining to the adaptive and innate immune responses.
Lastly, simulations are generated to make predictions, and statistical analyses are performed to extract useful information.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020. 10.30.20215335 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; bound at 10 2 copies ml -1 and 10 4.6 pg ml -1 is imposed to represent the detectable limit of viral load . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint and antibodies, respectively; dashed blue line indicates the limit of detection. Once the viral load or antibodies go below the detection limit, a vertical line is used to indicate the time of occurrence of this event. The triangular markers represent clinical data presented as mean ± SD (n = 4 patients). Pearson correlation coefficient between observed and fitted values is R > 0.98 (see Figure   S1b ).
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. Insets in a,d,g) show plasma concentration kinetics of the pharmacological agent/s for the treatment regimen that begins on the day of onset of symptoms. A lower bound at 10 2 copies ml -1 is imposed to represent the detectable limit of viral load in the body (black . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint horizontal line); once the viral load goes below the detection limit, a vertical line is used to indicate the time of occurrence of this event.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint Reequilibration rate of naïve CD8 cells cell⋅mL 1 ⋅d 1 -0.0045 62, 63 Reequilibration rate of naïve B cells cell⋅mL 1 ⋅d 1 -1.75 76 Reequilibration rate of naïve CD4 cells cell⋅mL 1 ⋅d 1 -0.0027 63, 62 Baseline IFN production rate IFN units⋅ mL -1 ⋅d -1 38.8 1.37 29, 30 Degradation rate of IFN d -1 24 24 77 Healthy tissue parameters Initial concentration of target cells in URT cells⋅mL -1 10 6.9 10 6.7 72, 78, 79 Initial concentration of target cells in LRT cells⋅mL -1 10 6.9 10 6.7
Initial concentration of target cells in GI cells⋅mL -1 10 7.4 10 6.9 80 Initial concentration of target cells in MPS cells⋅mL -1 10 6.8 10 6.6 81, 82, 83 Initial concentration of target cells in heart cells⋅mL -1 10 6.5 10 6.3
Initial concentration of target cells in kidneys cells⋅mL -1 10 7.1 10 6.9
Initial concentration of target cells in brain cells⋅mL -1 10 5.9 10 5.7
Viral dynamics parameters Initial concentration of viral load in URT virus⋅mL -1 10 5.7 Table 3 29
Clearance rate of antibodies d -1 -0.04 10 Hepatobiliary excretion rate d -1 2.68 0.15 89 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint Table 2 . List of model parameters estimated from fitting the reduced model to in vivo data.
Note: Characteristic times corresponding to the transport parameters are shown in Figure 1 . . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint Transition rate of naïve CD8 into CD8 * or CD4 into CD4 * mL⋅cell -1 ⋅d -1 0.0061 Transition rate of naïve B cells into B * mL⋅cell -1 ⋅d -1 0.029 Transition rate of B * into short-lived plasma cells mL⋅cell -1 ⋅d -1 564.8 Transition rate of B * into long-lived plasma cells mL⋅cell -1 ⋅d -1 0.36 Production rate of antibodies from short-lived plasma cells.
antibody⋅cell -1 ⋅d -1 0.59 Production rate of antibodies from long-lived plasma cells.
antibody⋅cell -1 ⋅d -1 1.34 Antibody loss rate due to viruses mL⋅virus -1 ⋅d -1 0.0013 Elimination rate of viruses due to antibodies mL⋅antibody -1 ⋅d -1 0.0085 Death rate of infected cells due to CD8 * mL⋅cell -1 ⋅d -1 0.0037 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)
The copyright holder for this preprint this version posted November 3, 2020. ; https://doi.org/10.1101/2020.10.30.20215335 doi: medRxiv preprint