COVID-19 is a recently emerging infectious disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [28] . After outbreaks of severe acute respiratory syndrome coronavirus (SARS-CoV) in 2002, and of Middle East respiratory syndrome coronavirus (MERS-CoV) in 2012, this is the third outbreak of a coronavirus since the turn of the century. Mathematical models of COVID-19 transmission at the population level have been instrumental in controlling the spread of the virus, but a detailed understanding of within-host infection dynamics is still lacking. Like SARS-CoV and MERS-CoV, SARS-CoV-2 is a betacoronavirus in Group IV of the Baltimore classification of viruses. Compared with other single-stranded (ss)RNA viruses, coronaviruses have the longest genomes. The smallest RNA viruses, for example, are only ∼ 1-4kb in length, and HIV has two copies of a ∼ 10kb genome, whilst the SARS-CoV-2 genome is a positive-sense, ssRNA molecule of ∼ 30kb. As a consequence, the viral life cycle of coronaviruses is distinct from most other ssRNA viruses, and existing intracellular infection models, e.g. for hepatitis C virus [2] , cannot be applied here. Therefore, we introduce here a novel intracellular model of SARS-CoV-2 infection, that incorporates essential steps specific to coronaviral life cycles. This model enables us to study in detail the viral dynamics inside an infected cell, and provides a framework to study the impacts of antiviral treatments on the infection dynamics within an infected cell.
In particular, it enables us to quantify the impact of different treatment options on the viral load that is secreted from an infected cell.
Outcomes from the intracellular model are then integrated into an intercellular model, that takes the impact of the immune response on infection dynamics within an infected individual into account. The model has been parameterised with data from 12 patients from a study in Singapore [41] , enabling us to generate a generic profile of disease progression in patients that have recovered from the disease. Model predictions agree well with experimentally and clinically measured parameters such as the duration of the incubation period, suggesting that this scenario is representative of disease progression seen in COVID-19 patients. We then use this model to study the infection dynamics in patients with different levels of immune responses by varying parameters associated with the immune response, such as the proliferation rates of effector cells and antibodies, and the rate by which effector cells remove infected cells. Comparison of different scenarios is based on tissue damage and viral load, highlighting the impact(s) of antibodies and adaptive cell-mediated immune response on infection dynamics.
This provides a framework in which to compare the impacts of different forms of antiviral therapy and assess their synergies. We focus here on two prominent forms of therapy against COVID-19: remdesivir, that inhibits virus production within an infected cell, and convalescent plasma (CP) therapy, whereby CP derived from recently recovered donors is transfused to the patients as an additional support [12] . As expected and previously observed [17, 25] , an early treatment start before the first peak in viral load can reduce both tissue damage and the peak viral load, especially when using a combination of both therapies. This reduction is quantified by the area under the viral load curve (AUC), which is an infection-related quantity commonly used in the assessment of a treatment against acute viral diseases [37] . Surprisingly, our model also suggests that early treatment by either therapy alone can actually increase the duration of infection compared with a later therapy start, likely because suppressing virus production results in a reduced immune response. This implies that even though early treatment accelerates the recovery process and reduces the peak in viral load, the infection may persist for a longer time than later treatment, meaning that these patients may possibly still be infectious. These effect is negated by using a combination of both therapies, but the synergistic benefit of combination therapies declines for later treatment starts. Our intercellular model provides insights into the effects of different dosages and treatment starts in terms of infection-related quantities, particularly for CP, finding the optimal dose and time for CP therapy is still debated [12] . Our model suggests that whilst a later treatment start does not result in a significant reduction of tissue damage, it can still reduce the AUC, thus potentially improving the course of infection and likelihood of recovery . The model also reveals that if treatment is started early, then antiviral therapy using remdesivir is more beneficial than CP therapy. By contrast, for a later therapy start CP therapy shows the better progression. This suggests that the time at which a patient was infected should be considered when deciding on the appropriate therapy.
Our stochastic model of viral infection dynamics within an infected host cell tracks the different viral and cellular components required for formation of progeny virus. These include the structural proteins that make up the virus: the envelope (E) protein, the membrane (M) protein, and the spike (S) protein, as well as the nucleocapsid (N) protein (Fig. 1a) . N protein forms a complex with the genomic RNA (gRNA), and thus aids its compaction for ease of packaging within the viral envelope. The S protein binds the ACE-2 receptor on human cells and is therefore essential for cell entry. The model also includes non-structural proteins that are important for the viral life cycle, such as the replicase-transcriptase complex (RTC), and keeps track of the numbers of gRNAs (and subgenomic sgRNAs) at different stages of the replication process. These include both the original plus-sense template of the RNA molecules, as well as their negative-sense variants that arise transiently during transcription. There are nine negative-sense subgenomic RNAs (-sgRNAs), corresponding to different gene products. In particular, there is an individual one for each of the structural proteins, allowing the virus to produce these components in the different quantities required for formation of viral particles [3] .
The reactions modelling viral replication within the host cell are described in detail in the Supplementary Information (SI). Here we provide a brief summary of the main reactions. The SARS-CoV-2 genome encodes two polyproteins, pp1a and pp1ab. The former is translated from the open reading frame ORF1a, and the latter from the overlapping reading frames 1a and 1b [15] via a -1 ribosomal frameshift during the translational elongation step. The relative frequencies of occurrence of the ribosomal frameshift is a key mechanism of self-regulation of protein expression of the virus and hence an important parameter in our model, which is captured by the frameshift probability q [32] .
Proteins cleaved from pp1a and pp1ab form the RTC, which is then used by the virus for genome replication and production of the nine -sgRNAs. -sgRNAs are produced through discontinuous transcription [11] , where elongation of nascent -gRNA continues until the first functional transcription-regulating sequence (TRS) is encountered. A fixed proportion of RTCs will disregard the TRS motif and continue to elongate the nascent strand, while the remainder will halt synthesis of the nascent minus strand and instead synthesise -sgRNAs [35] . The nine TRS motifs in the SARS-CoV-2 genome correspond to the nine sgRNAs produced, hence the choice to elongate or terminate synthesis occurs up to nine times during the elongation process [26] . The -gRNA and -sgRNAs produce positive-sense RNAs by recruiting RTC. The sgRNAs are translated to form proteins using cellular ribosomes. In the final step, a gRNA and the structural proteins (S, M, N, and E) form a new virion according to an assembly reaction that takes the stoichiometry of the different viral components into account [3] , and the virus particle is then released from the host cell.
Stochastic simulations of the reactions were implemented using a Gillespie algorithm [16] , and the number of particles released over the course of 100 hours was computed as the average over 200 stochastic simulations. Parameter values used (SI table S1) are predominantly based on experimentally available data; for parameters for which no data were available, we ensured that our main conclusions are robust against their variation. In particular, the release rate of virions, following a time lag between infection of a host cell and its first release of viral particles, is constant and virions are secreted linearly (SI Fig. S3 ).
We note, however, that the length of the time lag, and correspondingly the total number of virions released, are affected by some of these parameters and therefore warrant a more detailed investigation. For example, increasing the ribosomal protein production rate (SI Fig. S3a ) or the RTC nucleotide association rate (SI Fig. S3b ) decreases the time lag and therefore increases the number of virions released. By contrast, variation of the half-life of RTC or the formation rate of RTC from the constituent proteins does not have a significant effect on the time lag (SI Fig. S3c and d) . Figure 1b shows the time lag to the release of the first virion from an infected cell as a function of the ribosomal frame shifting probability q, and of the RTC elongation probability r.
An earlier treatment start, especially during the latent period, is more effective. (a) Illustration of a SARS-CoV-2 virion; the viral genome (dark green) is in complex with nucleocapsid (N) protein (light green) and is enclosed by the viral envelope that is studded by other structural glycoproteins, the spike (S) protein (red), the membrane (M) protein (maroon), and the envelope (E) protein (yellow). (b) Time lag before the release of the first virion from an infected cell; the maximal release of virions occurs when the RTC elongation probability, r, is high and the frameshifting rate, q, is between 0.3 and 0.5. (c) and (d) Profiles of viral load from an infected cell after introducing treatment at different times post infection, for a concentration of 25 (c) and 50 (d) molecules of remdesivir, respectively. The black curve indicates the drug-free control. The magenta, green, red, and orange curves correspond to a treatment start at 50, 30, 20, and 10 hours post infection, respectively. Parameter values are given in (SI Table S1 ).
This figure indicates that decreasing r results in a rapid increase in the time lag: for example, when r = 0.55 and the virus produces many sub-genomic fragments, this time lag is longer than 200 hours. This implies that viral load is maximised for the scenario that the RTC favours continuation of the transcription process when encountering a TRS. Given the importance of frameshifting in the coronavirus life cycle to control the relative numbers of different viral components, it has been argued that the virus will have evolved to optimise this ratio. In particular, frameshift signals have been characterized experimentally previously to have efficiencies in the range of 20-45% [4, 6, 22, 34] , although a recent study has suggested that the frameshift rate may be slightly higher [23] . Our model identifies an optimal value of 0.3 < q < 0.5 as the range with the lowest time lag and hence maximal virion release (Fig. 1b) , in good agreement with the experimentally determined range of 20-45%.
In order to assess the impact of an antiviral drug on viral load in the context of the intracellular model, we use the example of remdesivir, which is a widely used treatment option against COVID-19. Remdesivir, originally developed as a treatment for Hepatitis C virus and later trialled for efficacy against Ebola, acts as a nucleoside analogue that mimics adenosine nucleotide [38] . During the replication process, RTC may insert remdesivir molecules instead of adenosine, resulting in capping of the strand and thus terminating replication [42] . We have included additional reactions into the model that describe remdesivir binding to the RTC complexes on the gRNAs and sgRNAs to capture this (see SI for details), and track the effect of a given, fixed number of remdesivir molecules per cell on the release of viral particles from an infected host cell. Figures 1c and d show the impact on viral load released from a single infected cell for two different concentrations of remdesivir, as well as treatment starts at different times post-infection (TPIs). In Fig.  1c , a free concentration equivalent to 25 remdesivir molecules per cell is considered, corresponding to a concentration of ∼ 0.06 µM (see SI) [18] . The black curve indicates the viral load in the absence of treatment as a control, computed as an average over 200 stochastic simulations of the model. Magenta, green, red, and blue curves show the impact(s) of treatment start at TPIs of 50, 30, 20, and 10 hours, respectively. Figure 1c demonstrates that starting treatment during the latent period reduces the total number of virions released significantly. Even given a later treatment start the rate of virion production is slowed down, but the earlier treatment is started, the stronger the reduction in the virion production rate. Our results are consistent with experiments that probed the impact of remdesivir on MHV infection [1] , which also revealed that starting treatment earlier and during the latent period is more effective, as is the case also in other betacoronviruses. Figure 1d shows the impact of doubling the drug concentration (equivalent to a concentration of 50 molecules of remdesivir per cell). In this case, starting treatment early during infection at 20 hours post infection reduced the number of virions released on average by more than half. However, starting treatment later in the infection, such as 30 or 50 hours TPI, decreases the number of released virions by a smaller fraction. This suggests that although an increased drug concentration can be beneficial, starting the treatment earlier is more effective at reducing viral load than an increase in dosage.
The intracellular model affords insights into the release of viral particles from an infected cell, both in the absence and presence of antiviral treatment. We integrate results from this model into an intercellular model of within-host infection dynamics in order to probe the impact of the adaptive immune response on disease progression both in the absence and presence of antiviral therapy. Uninfected target cells (T ) are assumed to follow logistic growth with proliferation rate r T and carrying capacity T m . Inclusion of growth capacity of uninfected cells is important, because SARS-CoV-2 is detectable in patients over 20 days after the onset of symptoms [43] , comparable to the time taken to regenerate the epithelium (up to 1 month [39] ). Uninfected cells are infected by free virions at rate β . Although infected cells in the latent phase probably die at a somewhat lower rate than productively infectious cells, we assume that all infected cells die at approximately the same rate, δ , in order to minimise the number of free parameters that would complicate parameter estimation. Our intracellular model shows a time lag between infection of a host cell and the first release of new virions, consistent with experimental observation [3] . This effect is included into our intercellular model via a latent phase (L) with a lifetime defined as 1/γ, where γ denotes the average transition rate from the latent to the productively infectious (I) state, i.e. when the cell sheds viral particles. The intracellular model also shows that virions are effectively released linearly in time. Therefore we model infected cells as producing new virions V at a constant production rate p, and assume that they are naturally cleared at rate c.
Our model of the adaptive immune response consists of antibodies A (humoral immune response) that remove virions at rate k, and effector cells E (cell-mediated immune response) that kill infected cells at rate µ, assuming the same rate for cells in the latent and infectious phase in order to minimise the number of free parameters. Antibodies are produced at rate p A proportional to the viral load and are degraded at rate d A . After viral clearance, the antibody level is kept at a homeostatic level, because of the long-lived plasma and memory B cells. To represent this, we add a logistic term with proliferation rate r A and carrying capacity A m to the antibody equation. A fixed basal level of effector cells is assumed (λ E /d E ), and upon infection the population of effector cells will expand at rate α(L + I)E [7, 20, 14] . L and I both have an impact on the immune response, because infected cells during the latent phase are producing viral proteins. Although infected cells at different stages of infection are likely to express slightly different levels of viral peptide-MHC on their surface, we assumed that the rates are the same for L and I in order to minimise the number of free parameters in our model.
Considering the above assumptions, the model, as illustrated in Fig. 2 , takes on the following form:
(1)
The model was fitted to data from 12 hospitalised patients in Singapore [41] using measurements of viral load (see Materials and Methods). The parameter values derived from fitting V are presented in Table 1 . Our model captures essential features of the viral load in all patients, including the positions and heights of the first peak, and where applicable also those of the second peak (see SI Fig. S4 ). In all patients the viral load eventually decreases to below detectable levels, matching the clinical outcomes in these patients. We note that even details such as the slower viral decline in patients 2 and 12 are correctly represented by our model.
Our model predicts an incubation period, i.e. time between infection and presentation of symptoms, of 4.25 days (3.45-5.05 95% CI), in excellent agreement with the median SARS-CoV-2 incubation period of roughly 5 days estimated based on clinical data elsewhere [3, 29] . The average time after which antibodies appear is predicted here to be 16 days (13.9-18.1 95% CI) after infection, again in excellent agreement with the clinically reported first detection of antibodies after 10-20 days [3] . Similarly, the latent period of 27.28 hours (26.19-28.37 95% CI) predicted by our model agrees well with the experimentally observed latent period of 12-36 hours [21] .
The within-host model enables the roles of different aspects of the adaptive immune response in viral clearance to be investigated in more detail. The adaptive immune response to a viral infection relies on both antibodies and effector T cells. In order to understand their respective contributions to viral clearance, we first generate a generic progression profile based on the data from all 12 patients, and then vary parameters pertinent to different aspects of the immune response in isolation in order to probe their impact on disease progression.
The median values from our parameter fitting were used to generate a generic progression profile from the 12 patient data as a control (black curves in Fig. 3a ,b,c and SI Fig. S5 ). This control curve reveals a characteristic two-peak behaviour for viral load (Fig. 3b) , with antibodies passing the detection limit (0.1 ng/ml = 4 × 10 8 molecules/ml [8] ) after 14 days post infection (SI Fig. S5d) .
The cell-mediated immune response is captured in the equations by two factors: µ, the removal rate of infected cells by effector cells (T cells); and α, the proliferation rate of the effector cells. As λ E /d E is the basal level of effector cells, it is assumed to be constant for each patient [7] . Thus, α and µ are parameters that are likely varying in different patients. In particular, they would be expected to be lower in immunocompromised patients than for a patient with a healthy immune system [7, 19, 9, 14] . Figure  3d (red curve; see also SI Fig. S5 ) demonstrates that although reduction in the value of µ can increase the damage to healthy cells (T ) slightly, this effect is much stronger when reducing α (Fig. 3g , blue curves). Even though decreasing either of these parameters causes a slower decline in viral load after the second peak ( Fig. 3e and h) , a smaller value of α in addition increases the maximum of both peaks. Figures 3j,k,l (magenta curves in SI Fig. S5 ) model the case where just the humoral immune response is weakened, i.e. where the proliferation rate of antibodies r A is reduced [8] . In this case, viral load shows three peaks, and the damage to healthy cells recovers only slowly. This demonstrates that each component of the immune response plays a different, and crucial, role in the recovery process. In particular, in immunocompromised or elderly individuals who have a weakened immune response, this could lead to significant tissue damage, with infections lasting much longer than for non-immunocompromised patients.
In order to study the effects of antiviral therapy in the context of our intercellular model, we multiply the viral production rate p by (1 − ε), where 0 ≤ ε ≤ 1 is the drug efficacy [10, 27, 19, 14] . As our intracellular model ( Fig. 1c and d) suggests that starting treatment in the latent period is most effective, and would effectively block the production of virions, we set γ = 0 at the onset of treatment. This means that infected cells in the latent phase (L) do not transition to phase I, and begin shedding virions at a reduced rate (following a time lag τ) compared with cells that were already in phase I at the onset of treatment. We thus have the following equations for the numbers of infected cells and free virions; 
where θ (.) is the Heaviside function, τ denotes the time delay, and t R is the time when the antiviral treatment (remdesivir) is introduced.
As revealed by the intracellular model, the time at which the drug is administered impacts on the duration of the latent phase of an infected cell, which we model here using a time delay (τ). After this delay, infected cells produce virions at the reduced rate of p(1 − ε). Our intracellular model ( Fig. 1c and d) suggests that starting therapy during the latent period can reduce the level of virions released by ∼ 99% on average. We thus assume ε = 0.99 for the efficacy of remdesivir in cells in phase L. This is a good approximation, as all values of epsilon calculated are above 98%, regardless of treatment start time. Cells that were already in the productively infectious phase (I) at the time of treatment start are assumed to produce virions at rate p(1 − η), where η ≤ ε is the efficacy of remdesivir in this group of cells. Figure 1c indicates that starting therapy during the productively infectious period can reduce the level of released virions by ∼ 90% on average. Therefore we set η = 0.9. We note that these values of ε and η are consistent with values used in previous models [25, 26, 17] . As our numerical results show that the delayed model (SI Fig. S6 ) has the same behaviour as the model without a time delay, i.e. τ = 0, (Fig.  3) , we set τ = 0 in the following. Figures 3a,b ,c illustrate the impact of treatment on viral clearance in a patient who would most likely recover without treatment, as the parameters have been chosen as the median values of the 12 patients who have recovered from COVID-19 without treatment (Table 1 ). In this generic scenario, a treatment start 7 days post infection (dpi), which is approximately the time of the first peak in viral load, prevents the second peak from occurring (black dashed lines in Fig. 3 ). An earlier treatment start at 6 dpi also leads to the same result. However, starting treatment even earlier at 5 dpi (black dash-dot lines in Fig. 3 ) also reduces the damage to healthy cells. Even though free virus declines slower in this case, the area under the viral load curve (AUC), an infection-related quantity commonly used to help the assessment of a treatment against acute viral diseases [37] , is much smaller. In respiratory infections, even after viral clearance the immune response can cause respiratory and systemic symptoms in some incidences [37, 24] . A treatment start at 5 dpi results in a reduction in the peak of immune response cells, suggesting that early treatment perhaps could mitigate against this. Table 1 . Red curves correspond to the scenario of low removal rate of infected cells by effector cells (µ = 3.5 × 10 −4 ). Blue curves illustrate the scenario of a low proliferation rate of effector cells (α = 5.4 × 10 −10 ), and magenta curves of a low antibody proliferation rate (r A = 1). The green line indicates the viral detection limit.
In Fig. 3 solid red, blue and magenta curves indicate cases where different aspects of the immune response are weakened in isolation. In particular, the solid red curves illustrate the case of a reduction in the removal rate of infected cells by effector cells µ by 99% reduction with respect to the generic case (µ = 3.5 × 10 −4 ), the solid blue curves that of a 92% reduced proliferation rate of effector cells α (α = 5.4 × 10 −10 ), and the solid magenta curves correspond to a low antibody proliferation rate (r A = 1 instead of 1.98). As in the generic case above, figures 3d-l indicate that in each case starting treatment 5 dpi reduces tissue damage, viral peak height and AUC significantly, compared with treatment starts at 6 or 7 dpi, again emphasizing the importance of an early treatment start. However, early treatment increases the duration of infection compared with a later therapy start (Fig. 3h) . This suggests that viral load could persist for a longer time in such patients, who may still be infectious.
Our intercellular model also enables the modelling of drugs that operate at the level of the immune response, rather than virus production in the intracellular milieu. As an example of a therapy option of this type we study the impact of convalescent plasma (CP) therapy on viral dynamics [12] . For this, we add a new variable to our model (Ã), which captures the antibodies that are administered as treatment. We assume that these antibodies remove virions at rate k f , where 0 ≤ f ≤ 1, implying that they are at most as efficient as the antibodies that are being developed by the body over the course of infection. The equations for the number of virions and antibodies thus have the form;
whereÃ(t) = 0 for t < t CP andÃ(t CP ) =Ã m , withÃ m representing the number of antibodies per ml that are administered as treatment. t CP denotes the time at with the treatment is started, and θ (.) is the Heaviside function.
Our model enables the impact of CP therapy on viral dynamics to be studied for different treatment starts and doses, thus addressing the bottle-neck pointed out in the recent literature of finding the optimal dose and start for CP treatment [12] . Using again median values from Table 1 to generate a generic patient profile as a control, and using three immunocompromised cases that are presented in Fig. 3 (red, blue, and magenta curves, representing cases with reduced values for the immune response parameters µ, α, and r A , respectively) we studied the impact of CP therapy. SI Fig. S7 shows the minimum level of therapeutic antibodiesÃ m that is needed to reduce the AUC by 25% and 50% as a function of the start of treatment (in dpi) and the factor f by which therapeutic antibodies are less efficient than those produced by the host during the infection. It indicates that the level ofÃ m that is needed for an effective reduction in the AUC is at most around 3 × 10 11 molecules/ml. This is in good agreement with clinical data, reporting a 200 ml dose of CP and A m = 4 × 10 12 molecules/ml (see Materials and Methods) [12, 36, 8] . Indeed, this implies that each dose would contain about 8 × 10 14 molecules, resulting iñ A m = 2.6 × 10 11 molecules/ml on the basis of an average level of 3 liters blood in the body [31] . Hence, a reduction of the AUC by 50% is achievable. SI Fig. S7 also shows that AUC reduction is comparable in the range 0.7 ≤ f ≤ 0.9, therefore we use the average value of this range ( f = 0.8) in our calculations. Our conclusions are robust for efficiencies ≥ 0.15, while for values of f below 0.15, the outcomes vary for different immunocompromised cases. Using these parameters, we present a comparative analysis between antiviral and CP therapy and explore their synergistic potential. Figure 4 shows that similar to the case of antiviral therapy in Fig. 3 , an early treatment start is more effective in the reduction of tissue damage and the level of AUC compared with a later therapy start (cf. SI Fig. S8 for the equivalent of Fig. 3 for CP therapy). While Fig. 3 reveals scenarios in which an antiviral therapy does not mitigate against tissue damage (such as a later treatment start at 6 or 7 dpi), Fig. 4 shows that using CP therapy can reduce tissue damage even for those delayed treatment starts. Interestingly, Fig. 4 shows that starting CP therapy early can increase the duration of infection more than for antiviral therapy, implying that for an early treatment starts using remdesivir is more effective. By contrast, for later treatment starts CP therapy reduces the viral load faster and decreases tissue damage compared with remdesivir therapy. Our model also enables us to probe the synergies of these treatments options. Figs. 4c, f, i, and l indicate that for an early treatment combination therapy mitigates against a longer duration of the infection. However, for later treatment starts any synergistic effects are minimal and combination therapy has the same outcome as CP therapy in isolation. Thus, unless infection is detected early, e.g. though an efficient track and trace system, treatment would likely start at a time when CP therapy in isolation would be as effective as combination therapy. Table 1 . Red curves correspond to the scenario of low removal rate of infected cells by effector cells (µ = 3.5 × 10 −4 ). Blue curves illustrate the scenario of a low proliferation rate of effector cells (α = 5.4 × 10 −10 ), and magenta curves of a low antibody proliferation rate (r A = 1). The green line indicates the viral detection limit. First and second columns indicate the impact of only CP therapy while third column show the impact of having both antiviral and CP therapy. SI Fig. S9 indicates the impact of starting treatments after the onset of symptoms corresponding to the 12 patients used for model fitting, both treatments would have reduced the duration of the infection significantly (in 67% of the patients), enabling a faster recovery, or otherwise has no significant impact on the duration of the infection. However, in some cases although the peak in viral load, the AUC and tissue damage are significantly reduced, the treatments have not decreased the duration of infection. This figure also shows for the cases where the duration of infection is not reduced by one treatment, using both treatments would have reduced it. In the other cases these treatment options have more or less similar effects although CP therapy acts slightly better and there is not a noticeable synergistic effect between them.
The severe consequences of the COVID-19 pandemic demand a concerted interdisciplinary effort to identify novel antiviral solutions. Vaccination options against SARS-CoV-2 are actively pursued and a number of treatments are now in use in the clinics, but there are still many open questions regarding when and how best to administer these treatments, either in isolation or in combination. Whilst modelling of disease transmission has already played a key role in informing policy makers [5] , models of within-host dynamics have not yet had a prominent role in combating the disease. There is a precedence of intracellular modelling for other viral diseases, such as hepatitis C virus [2] . However, such models cannot be readily transferred to coronaviral infection, as viral life cycles are very different. Here we introduce a within-host model of a SARS-CoV-2 infection that contains sufficient details specific to coronaviruses to enable antiviral strategies against SARS-CoV-2 to be compared and to analyse their synergies. We demonstrate this via a comparative analysis of an antiviral treatment (remdesivir) and CP therapy, which apart from steriod treatment are the most prominent forms of therapy currently used against COVID-19 infections. In particular, we compare disease progression for different treatment starts and dosages, and thus provide new insights into these therapeutic options.
As expected, our analysis highlights the importance of early treatment starts post infection, consistent with other model predictions [25, 17] . However, those models do not capture the impact of early treatment on the duration of the infection. Our model has provided insights into disease progression for different doses and treatment starts for CP therapy option [12] . In particular, it enabled us to address a question recently raised in the literature as to the impact of dose and time of treatment on disease progression under CP treatment [12] . Our model also enabled us to perform a comparison between the antiviral treatment and CP therapy, and explore their potential synergistic effects. The model reveals that early into the infection an antiviral treatment using remdesivir could be more effective than CP therapy, and a combination therapy can significantly reduce the duration of infection. However, later for later treatment starts, CP therapy appears to be more beneficial than antiviral therapy, and there are no longer any significant synergistic effects that would warrant combination therapy. These insights from our within-host model suggest that the time course of infection should be considered when deciding on appropriate theraputic response to COVID-19 infection.
Our patient comprise the first 18 confirmed patients who reported COVID-19 infection in Singapore [41] . Nasopharyngeal swabs were collected for up to 30 days since onset of symptoms. Five patients received lopinavir-ritonavir treatment, and in one patient viral load was detectable only twice, and these six patients were therefore excluded from the analysis. The viral loads were reported in cycle threshold (Ct) values, which is inversely proportional to viral RNA copy number [43] . We converted Ct values to viral copies per ml. In model fitting, viral load values under the detection threshold were set at the detection limit (Ct=38).
COVID-19 is a respiratory illness, so we assume that modelling insights from influenza models are applicable. In influenza, at approximately 5 to 7 dpi mitoses are detected at the basal cell layer, and regeneration of the epithelium begins. Complete resolution of the epithelial takes up to 1 month [39] . We therefore assume that the maximum proliferation rate for uninfected cells is small, and that r T = 0.1 day −1 . The number of host cells that express ACE2 and TMPRSS is approximately equal to 10 11 (T m = 10 11 ) [3] , and we use T (0) = T m . As SARS-CoV-2 is a novel infection, we assume that E(0) = 0 and that the basal level of effector cells is low λ E = 1 and d E = 0.5 day −1 [8] . However, considering a higher basal level does not change model outcome as α and µ which determine the proliferation rate of effector cells and removal rate of infected cells by effector cells are estimated using data fitting. Initially there is no specific antibody, therefore A(0) = 0 and d A = 0.033 day −1 [8] . We use 1 µg/ml IgG positive control as a strong positive standard [36] . Thus, we assume A m = 1 µg/ml = 4 × 10 12 molecules/ml [8] .
The patient data used is only available from the time after onset of symptoms, and the initial viral load at the start of the infection is not recorded. We therefore estimate the value V (0) assuming the infection is transmitted via droplets. The average number of expelled droplets during talking is assumed to be 1000 [30, 40] . It has also been reported that more than 50% of droplets have a size range between 50-75 µm [40] . Thus, the average volume in expelled droplets during talking is equal to 1.1 × 10 −4 ml. The median level of viral load on the day of symptom onset in patients in this study is estimated as 5 × 10 3 virion/ml [13] . We assume that infected individuals infect others before the onset of symptoms, and we therefore assume an average level of 10 3 virion/ml are available for transmission. Thus, assuming V (0) = 0.1 virion/ml appears to be a reasonable choice [17] . This value is also comparable to those used in modelling of influenza [33] .
We estimate the remaining parameters and the incubation period (the time between the beginning of the infection and the onset of symptoms) by fitting V from model (1) to patient data using MONOLIX 2019R2 (www.lixoft.com), a program for maximum likelihood estimation in nonlinear mixed-effects models. The resulting parameter values are presented in Table 1 . Decreasing/increasing of V (0) (V (0) = 0.01 virion/ml or V (0) = 1 virion/ml) does not change the estimated values of parameters significantly and only changes the estimated incubation period by ±1 day.
Farzad Fatehi, Richard J Bingham, Eric C Dykeman, Peter G Stockley, and Reidun Twarock * * Corresponding Author: rt507@york.ac.uk S1 A stochastic model of intracellular dynamics of SARS-CoV-2
The intracellular model contains reactions representing the translation, transcription and assembly of SARS-CoV-2 in cells. The viral entry step is not included into the model as we are focusing on the intracellular dynamics. Here we will detail the reactions underlying each process and describe the stages which have been modelled mathematically. All parameter values are given in Table S1 . Stochastic simulations of the reactions were implemented using a Gillespie algorithm [1] , and seeded with one virion at the start of the infection.
The SARS-CoV2 genome, illustrated in Fig. S1a , shows the gene products encoded by the 29,903 nts ssRNA viral genome. The replicase gene comprises roughly 22kb and encodes two polyproteins (pp1a and pp1ab) which undergo cotranslational proteolysis into the 16 non-structural proteins (nsp1-nsp16) required by the virus for infection in the host cell. The remaining genes 3' of the replicase gene encode for the structural proteins S, M, E, N. Structural and non-structural proteins are translated from different messenger RNAs that are produced by the viral replicase-transcriptase complex (RTC) from viral gRNA. The replicase gene is translated from full-length viral RNA whereas the structural genes are translated from sub-genomic mRNAs synthesised by the viral RTC [2] . We model the translational process from the mRNAs produced by the RTC as follows (Fig. S1b) . For synthesis of the polyproteins pp1a and pp1ab, host ribosomes (denoted by R in Fig. S1b ) reversibly bind to with binding/unbinding rates R on and R of f , K d = 0.08 µM [3] . As the ribosome can bind and unbind, we assume after binding it can either unbind or start the elongation process. Thus, we model the kinetic steps involved in ribosome initiation and transition to the elongation state (Ri 1a :(+)RNA) that occur subsequent to ribosomal binding to the as a single kinetic step with rate R in . The ribosome then translates the pp1a gene (ORF1a) at rate T 1a . After translation of the pp1a gene, the ribosome can either frameshift -1 nt to the ORF1b reading frame, translating the polyprotein pp1ab, or terminate, releasing the polyprotein pp1a [4] . We model the -1 ribosomal frameshift as a reaction with rate qT f and the termination at ORF1a as a reaction with rate (1 − q)T f . These rates give a probability of frameshifting from the ORF1a to ORF1b reading frames of q. If the ribosome sucessfully frameshifts, it completes translation of the ORF1b reading frame with rate T 1b , terminates, and releases the polyprotein pp1ab. As the viral genome needs to bind to the RTC for replication, we assume that (+)RNA can have only one ribosome or one RTC on at a time. Details of the reactions are listed below. The polyproteins pp1a and pp1ab are cleaved by two proteases, papainlike protease (PLpro; corresponding to nsp3) and a main protease, 3C-like protease (3CLpro; corresponding to nsp5). Polyproteins pp1a and contains the nsps 1-11 and pp1ab contains nsps 1-10 and nsps 12-16 [2, 5] . It has been proposed that nsp2-16 collectively constitute a functional RTC in infected cells [6, 7] . These nsps have their own specific roles in the replication process, but the functions of some of them are not well understood [8] . For example nsp12 is the RNA-dependent RNA polymerase (RdRp); nsp13 is the NTPase/helicase; nsp14 is a proof-reading exonuclease [2] . As the exact stoichiometry of these nsps in an RTC molecule is unknown [9] and nsp11 which is a part of RTC is only in pp1a, we assume that pp1a and pp1ab have equal impact on RTC formation and this step is modelled as follows:
pp1a + pp1ab [10] . The other (+)sgRNA encode accessory proteins which are hypothesised to interfere with the host innate immune response and their functions are poorly understood [11] . Thus we only include Synthesis of S, M, N, and E proteins into our model. Furthermore, similar to previous work on intracellular modelling of HCV, where the number of free ribosomes is fitted to obtain the viral dynamics that was observed experimentally [12] , we assume that free ribosomes are at an equilibrium level, where sgRNAs are saturated with available ribosomes and produce protein at a constant rate T n , T m etc. The E protein is found in small quantities within the virion. We assume each SARS-CoV-2 virion has 100 copies of E protein assembled into 20 pentameric structures, and 300 copies of S protein assembled into 100 trimeric structures [13] . The level of M and N proteins in SARS-CoV-2 consists of approximately 2000 and 1000 copies, respectively [13] . Note that in this model we assume that SARS-CoV-2 has the similar levels of E, S, M and N proteins to SARS-CoV-1 which was reported for [13] . Similar to the previous intracellular models for HCV infection which have modelled the assembly and budding of virions as a single reaction [12] , we model the budding of a virion is modelled as a single reaction with budding rate k bud as follows:
Transcription of the RNA genome of SARS-CoV-2 occurs via the interaction of the RTC complex with genomic (+)ssRNA, genomic (-)ssRNA, and negative sense sub-genomic fragments. During transcription of (+)ssRNA, a subset of 9 sub-genomic minus strands, which encode all structural proteins, are produced through discontinuous transcription [11] . Each sub-genome contains a 5 leader sequence corresponding to the 5 end of the genome which allows interaction between host ribosomes and (+)sgRNAs. Each sub-genomic RNA consists of a single ORF that encodes for a structural or accessory protein [14] . In SARS-CoV-2, the order of these ORFs (from 5' to 3') is S, 3a, E, M, 6, 7a, 7b, 8, and N [10] . The full length (+)ssRNA genome contains functional transcription-regulating sequence (TRS) motifs which are found at the 3 end of the leader (leader TRS) and in front of each of the 9 ORFs. During transcription of the full length minus strand by the RTC, the process can terminate at one of these TRS motifs, resulting in one of the 9 negative sgRNA being produced. In our model, when an RTC encounters TRS motif number k, it will continue the elongation of the negative strand with rate rT c , and terminate with rate (1−r)T c , resulting in the production of (-)sgRNA k after extension by transcription of the 5 end of the genome [14] . The completed minus-strand sgRNA serves as a template for (+)sgRNA synthesis. Figure S2 illustrates the transcription reactions that we model in this work. For transcription from full length genome (+)ssRNA, RTC (denoted by RTC in Fig.  S2 ) reversibly binds to 3' UTR with binding/unbinding rates RT on and RT of f , K d = 0.1µM [15] .
We model the kinetic steps involved in RTC initiation and transition to the elongation state that occur subsequent to RTC binding to the 3' UTR as a single kinetic step with rate RT in . The formation of (-)sgRNAs can be modelled as follows: Here, T r k is the rate of RTC transcription between the TRS at site k − 1 and site k (Fig. S2) and each sub-genomic (-)RNA for k = 1, 2, ..., 9, corresponds to sgRNAs for N, 8, 7b, 7a, 6, M, E, 3a, S (-)sgRNAs, respectively, whereas the k = 10 state with transcription rate T 10 denotes the rate to transcribe the remaining length of 22kb of RNA upstream of the structural genes. This last step is responsible for the creation of full length (-)genomic RNA.
To model the synthesis of the 9 (+)sgRNAs that are required for translation by the host ribosomes, we assume that RTCs can bind to negative sense RNAs producing a polysome arrangement with up to n k RTCs present on (-)sgRNA k . RTCs bind to (-)sgRNA one at a time with binding rate b k , until the (-)sgRNA k is saturated with n k RTCs. After saturation, the last RTC (i.e. the n k th) at the 5' end of the (-)sgRNA can terminate transcription with rate p i , resulting in the production of a (+)sgRNA and release of an RTC. The resulting (-)sgRNA with n k − 1 RTCs can bind another RTC at the 3' end. This reaction models the movement of the remaining RTCs down the RNA strand and movement of the 5' most RTC to the 5' end for termination and daughter strand release. The synthesis of (+)sgRNA is described via the following reactions: 
Remdesivir acts as a nucleotide analogue that mimics the adenosine structure. It was originally developed as a treatment for Hepatitis C virus and later repurposed for Ebola [16] . Recent studies have pointed to remdesivir as an effective antiviral treatment option for Covid-19 [17] . During the replication process RTC may insert remdesivir molecules rather than adenine, which caps the strand and stops the replication process [18] . In order to model the impact of this drug, we assume that complexes with RTC in our model can bind (and subsequently unbind) to remdesivir molecules (Rem). Thus the new reactions have the following form: S1.4 Parameters for the intracellular model S1.4.1 RNA and protein decay rates Wada et al. [19] have measured half-lives of the genome and S-sgRNA for mouse hepatitis virus (MHV), a prototypic member of the CoV family belonging to the same genus (Betacoronavirus) as SARS-CoV-2. Half-lives of S-sgRNA and genome are measured as 5.96 and 5.41 hours, respectively. Therefore, we assume that half-lives of (+)RNA, (-)RNA, and all positive and negative sgRNA are 6 hours. Thus, their decay rate is equal to ln 2/6 hour −1 . Half-life of E protein is 1 hour [20] . We are assuming that half-lives of N, M, S, pp1a and pp1ab are also 1 hour. As half-life of the viral polymerase of HCV (like SARS-CoV-2 a +ssRNA virus) has been measured as 6 hours [12] , we assume that half-life of RTC is also 6 hours. See Table S1 for full details on these constants.
Dimelow and Wilkinson [3] have built a detailed kinetic model of translation initiation in yeast. They have estimated that the forward rate constants for ribosome binding are in the range of 1-50 (µM s) −1 . In this model we assume that the binding rate of ribosome to RNA is 25 (µM s) −1 and its unbinding rate is 2 s −1 [3, 21] . RTC binds to RNAs with dissociation constant (K d = R unbind R bind ) equal to 0.1 µM [15] . It has been observed that increasing/decreasing RTC binding rate does not have a significant effect on the outcome as it will also increase/decrease RTC unbinding rate [22] . Therefore we consider that RTC binding rate is equal to 1 (µM s) −1 and the unbinding rate is 0.1 s −1 . The translation rates are computed using the length of each protein and considering a constant speed for the ribosome. We study the impact of changing of this parameter (ribosome speed) in Fig. S3a . We observe that by changing ribosome speed the model behaves qualitatively the same. Thus, we fix this parameter in other simulations and we assume that the speed of the ribosome is 10 amino acids per second. The replication and transcription rates of RTC are computed using the length of each RNA and considering a constant speed for RTC [10] . We assume that the speed of the RTC is 30 nucleotides per second [23, 24] , although this parameter has also been varied. See Table S1 for full details on these constants.
The relative free energy of binding for remdesivir is ∆G = −8.28±0.65 kcal/mol [18] . The following equation relates the binding free energy (∆G) to K d :
, k B is the Boltzmann constant and T is Kelvin temperature [25] . Therefore, for remdesivir K d is approximately 1 µM and we assume that the binding rate is equal to 1 (µM s) −1 and the unbinding rate is 1 s −1 .
The binding rates are measured in (M s) −1 , but as we are using a reaction based model with the Gillespie algorithm, we use the cellular volume to change these units to (molecule hour) −1 [25] . For changing the units of parameters related to the ribosome, we have used the volume from yeast, but for RTC and remdesivir, we have used the volume of an E. coli cell as these parameters were measured in E. coli. The level of drug in each cell is constant and equal to 25 molecules. Since V = 0.7µm 3 , it is equivalent to 0.06 µM [26] . The ribosome translates the pp1a gene (ORF1a). After translation of the pp1a gene, the ribosome can either frameshift -1 nt to the ORF1b reading frame, translating the polyprotein pp1ab, or terminate releasing the polyprotein pp1a. Probability of frameshifting from the ORF1a to ORF1b reading frames is q. Figure S2 : Transcription of SARS-CoV-2 genome via the RTC complex. During transcription process a subset of 9 sub-genomic minus strands are produced through discontinuous transcription. The full length genome contains functional (TRS) motifs which are found at the 3 end of the leader (leader TRS) and in front of each of the 9 ORFs (Fig. S1a) . When an RTC encounters TRS motif, it will either continue the elongation of the negative strand or terminate the transcription, resulting in the production of (-)sgRNAs. Probability of RTC elongation is r. Table S1 . Total number of secreted virions from an infected cell as a function of (a) ribosome amino acid association rate, and (b) RTC nucleotide association rate. Total number of released viral particles as a function of (c) RTC half-life, and (d) RTC formation rate (f ). Table 1 (black lines are given by µ = 5.7 × 10 −2 , α = 7.08 × 10 −9 and r A = 1.98). A reduction in the removal rate of infected cells by effector cells causes a spike in the level of effector cells (red lines given by µ = 3.5 × 10 −4 ). A decrease of the proliferation rate of effector cells causes significant damage to healthy cells and increases the viral peak (blue lines are given by α = 5.4 × 10 −10 ), and reduction of the proliferation rate of antibodies indicates a third peak in viral load (magenta lines are given by r A = 1). Table 1 with τ = 2 days. Red curves correspond to the scenario of low removal rate of infected cells by effector cells (µ = 3.5 × 10 −4 ). Blue curves illustrate the scenario of a low proliferation rate of effector cells (α = 5.4 × 10 −10 ), and magenta curves of a low antibody proliferation rate (r A = 1). The green line indicates the viral detection limit.
Outcomes are similar to the model in the absence of time lag, which is therefore used in the main text. Figure S7 : The impact of CP therapy on the reduction of AUC for different treatment starts and antibody efficacy. The first and second column illustrate the minimum level ofÃ m required for 25% and 50% reduction in the AUC, respectively. (a) and (b) Parameters are the median values of Table 1 , i.e. represent the generic case based on the patient data; (c) and (d) correspond to the scenario of low removal rate of infected cells by effector cells (µ = 3.5 × 10 −4 ); (e) and (f) illustrate the scenario of a low proliferation rate of effector cells (α = 5.4 × 10 −10 ); and (g) and (h) show a case with a low antibody proliferation rate (r A = 1). Table 1 . Red curves correspond to the scenario of low removal rate of infected cells by effector cells (µ = 3.5 × 10 −4 ). Blue curves illustrate the scenario of a low proliferation rate of effector cells (α = 5.4 × 10 −10 ), and magenta curves of a low antibody proliferation rate (r A = 1). The green line indicates the viral detection limit. Figure S9 : Starting treatment after the onset of symptoms reduces the peak in viral load and leads to faster viral clearance. Solid lines indicate the best fit to patient data. Dashed and dotted curves indicate the result of starting antiviral and CP therapy after the onset of symptoms, respectively, while dash-dotted curves indicate a combination therapy start after the onset of symptoms. 