Success of prophylactic antiviral therapy for SARS-CoV-2: Predicted critical efficacies and impact of different drug-specific mechanisms of action

Repurposed drugs that are safe and immediately available constitute a first line of defense against new viral infections. Despite limited antiviral activity against SARS-CoV-2, several drugs are being tested as medication or as prophylaxis to prevent infection. Using a stochastic model of early phase infection, we evaluate the success of prophylactic treatment with different drug types to prevent viral infection. We find that there exists a critical efficacy that a treatment must reach in order to block viral establishment. Treatment by a combination of drugs reduces the critical efficacy, most effectively by the combination of a drug blocking viral entry into cells and a drug increasing viral clearance. Below the critical efficacy, the risk of infection can nonetheless be reduced. Drugs blocking viral entry into cells or enhancing viral clearance reduce the risk of infection more than drugs that reduce viral production in infected cells. The larger the initial inoculum of infectious virus, the less likely is prevention of an infection. In our model, we find that as long as the viral inoculum is smaller than 10 infectious virus particles, viral infection can be prevented almost certainly with drugs of 90% efficacy (or more). Even when a viral infection cannot be prevented, antivirals delay the time to detectable viral loads. The largest delay of viral infection is achieved by drugs reducing viral production in infected cells. A delay of virus infection flattens the within-host viral dynamic curve, possibly reducing transmission and symptom severity. Thus, antiviral prophylaxis, even with reduced efficacy, could be efficiently used to prevent or alleviate infection in people at high risk.


A. Continuous virus production model
We recall the model from the main text. Infectious and non-infectious virus particles are denoted by V I and V N I , respectively, target cells by T , infected cells in the eclipse phase by I 1 and infected cells producing the virus by I 2 . We use a previously studied within-host model of virus production (Pearson et al., 2011;Conway et al., 2013). The underlying individual-based reactions are the following: infection of target cells, Since we model the early state within-host dynamics of a viral infection, we can assume that the number of infectious virus particles, V I , is low so that the number of target cells is not strongly affected by transformation to infected cells, i.e. T (t ) ≈ T (0) = T 0 . Then, the first reaction can be rewritten as Using standard techniques to derive a set of ordinary differential equations from this set of reactions (e.g. Anderson and Kurtz, 2011), we find the system given in Eq. (1) in the main text. Note, that in the main text we use capital letters to denote densities while here the capital letters refer to the actual numbers of cells and virus particles.

A. . Connection to a burst model
Since the individual-based model is built on stochastic interactions of cells and virions, the number of virions produced by an infected cell is a random variable. Assuming that all virions are released at a single time, typically at cell death, the number of released virions, the burst size, follows a geometric distribution (Hataye et al., 2019). This can be seen by the following reasoning: the life-time of an infected cell is exponentially distributed with mean 1/δ and during this time there is a continuous production of virions at rate p. This production, assuming that it is a Markovian process, is described by a Poisson process (see Anderson and Kurtz (2011) for the general theory of modeling chemical reactions). The probability of the burst size, denoted by N , to be of size j is then given by the following calculation: This is the distribution of a geometrically distributed random variable with success probability p/(p + δ). Intuitively, the infected cell has undergone j + 1 steps, where the initial j steps resulted in the production of a virus (term (p/(p + δ)) j ) and the ( j + 1)-th step was its death (term δ/(p + δ)). The mean of this geometric distribution is p/δ. The continuous-production model can therefore be seen as equivalent to a burst size model with a burst size N having a geometric distribution with mean p/δ.

B. Burst model
The continuous-production model is more likely to be relevant for SARS-CoV-2 (Park et al., 2020), and was therefore chosen in the main text. Here we examine how a burst model would affect our findings. In a burst model, we assume that virus is produced in an infected cell but is only released to the environment upon cell death. The number of virus particles released is therefore a random number which we again denote by N .
In the corresponding reactions in Eq. (A.1), we need to replace the virus production and cell death lines by In order to be consistent with the continuous-production model, we set the mean of the burst size to p/δ. In the following we assume that the overall burst size N is Poisson-distributed. There are two reasons for this choice: (i) it is analytically relatively easy to handle, and (ii) it represents the other end of the spectrum of negative-binomially distributed burst sizes when compared to the continuous-production model which is equivalent to a geometrically-distributed burst size. The negative binomial distribution models the number of failures before a certain number r of successes in a sequence of Bernoulli experiments with success probability q. Specifically, the probability of observing k failures before the r th successes is given by The parameter r is sometimes also referred to as the dispersion and formally the negative binomial distribution can also be extended to non-integer values of r . The mean is given by qr /(1 − q). The negative binomial distribution relates to the geometric distribution by setting r = 1 and to the Poisson distribution by letting r → ∞. The probabilities of establishment for the continuous-production model and the Poisson-distributed burst size model therefore represent the two extremes of negative-binomially distributed burst size models with dispersion parameter r ∈ (1, ∞). This holds because the establishment probability can be computed by the probability generating function (Haccou et al., 2005) which is continuous and monotone in the dispersion parameter r . It is given by where z is an auxiliary variable.

B. . Establishment probability
We compute the establishment probability of the virus in the burst size model. A key ingredient is the offspring distribution of a single virus particle. The offspring distribution is given by a The life cycle of a virus (conditioned on infecting a cell) is given by a three step process: cell infection, eclipse phase and virus production within an infected cell. Ignoring this time delay which is irrelevant if we just consider the establishment probability, the virus population can be modeled by a discrete time branching process. At each time, all infectious virions alive at the time step before produce a random number of (infectious) virions according to the offspring distribution given in Eq. (B.4).
The extinction probability of a time-discrete branching process, when starting with one infectious virus particle, is given by the non-trivial fixed point of the probability generating function of the offspring distribution, i.e. the fixed point in the interval (0, 1) (Haccou et al., 2005). The probability generating function is given by where z is an auxiliary variable and E denotes the expectation of the random burst size of infectious virions ηN . The fixed point of this function is given as where W (x) is the Lambert-function (sometimes also called the product logarithm). It is defined for x ≥ − exp(−1). For values below this threshold, we need to solve Eq. (B.5) numerically. In fact, when plotting the establishment probability in Fig A below, we solve Eq. (B.5) numerically because the approximation of the Lambert-W function W (x) is inaccurate for negative x, especially when close to − exp(−1). The establishment probability, denoted ϕ, is then given by where V I (0) is the initial number of infectious virions. For alternative derivations of this result see also Pearson et al. (2011) and Conway et al. (2013).

C. Comparison of the continuous-production and burst model
We compare the establishment probability from the burst model described above with that obtained in the continuous-production model. Redrawing the first row of Fig 1 from the main text and comparing it with the corresponding graphs obtained from the burst model, we do not see any qualitative difference between the two models, cf. Fig A. As outlined in Section B, the two studied models can be seen as the extreme values of a model continuum.
By varying the dispersal parameter r of the negative binomial distribution, one can explore the entire continuum between the geometrically distributed burst size (which is equivalent to the continuous-production model) and the Poisson-distributed burst size model. Therefore, it seems safe to say that the exact mechanism by which we implement virus production in the model will only result in minor quantitative differences on the probability of virus establishment.

Parameter set
Low burst size (LowN) 11.2 4 × 10 4 18.8 7.69 High burst size (HighN) 112 4 × 10 3 188 7.69 Table A: Model parameters used in the main text and for the simulations in Fig A. The remaining parameters are not changed between the simulations and are set to:

D. Establishment probability when starting with a single infected cell
In this section, we investigate how the establishment probability changes if treatment is started when there is already an infected cell within the host. This situation might be more realistic to post-exposure treatment where infectious virus from the initial inoculum might have already infected a target cell (if the virus was not cleared). Instead of starting with a viral inoculum, we thus need to consider the situation where an infectious cell is already producing virus (but has not yet produced an infectious virus particle). The reasoning for computing the establishment probability is then as follows: we combine the establishment probability with initially j infectious virus particles with the probability for this infected cell to produce j infectious virus particles. As we have seen in Section A.1 the number of infectious virus particles produced by an infectious cell is geometrically distributed with success parameter δ/(δ + ηp). Therefore, the establishment probability when starting with an infected cell, denoted by ϕ I , is given by This result has also been derived in Pearson et al. (2011), where this analysis was done for the continuous-output and the burst model, and in Duwal et al. (2019) for a similar model in the context of HIV prophylaxis.
In our high burst size parameter set, there is no visible difference between treatment with a drug reducing productivity p and a drug reducing viral infectivity β (Fig Bb). However, for the low burst size parameter set, in contrast to what we found in the main text when initializing the system with a viral inoculum, now drugs reducing the infectivity p (blue) stronger reduce the establishment probability than drugs reducing the infectivity β (orange), cf. Fig Ba. This is explained by the order in which the drugs act: while a drug reducing viral production can immediately lower the chances for a further virus propagation, drugs reducing infectivity need to 'wait' for their targets, the extra-cellular virus, to arrive.

E. Combination therapy in the HighN parameter set
We investigate the effect of combination therapy in the high burst size parameter set (Table A). We find that the overall shape of the curves do not change compared to the LowN parameter set. A higher burst size decreases the establishment probability of the virus. If we compare Fig Cb with Fig 3 in the main text, we see that a ten-fold increase of the initial inoculum in the HighN parameter set (V 0 = 10) gives similar quantitative results as the LowN parameter set with V 0 = 1. This can be attributed to our ten-fold increase of the burst size when deriving the HighN parameter set from the LowN parameter set.

F. Time to detectable viral load
In this section, we study the mean time to reach a certain amount of viral load at the infection site within the host. We approximate this time using a mixture of deterministic and stochastic arguments. Classical branching processes typically have two possible outcomes: either the process goes extinct or grows indefinitely (Haccou et al., 2005). The deterministic model is captured by the mean of such a branching process, i.e. it takes into account both possible outcomes. Therefore, if we condition the branching process on survival, the deterministic model will typically underestimate the actual size of the corresponding branching process (Desai and Fisher, 2007). One can correct this error by rescaling the deterministic process by the probability of survival. In our specific setting this means that the total number of virus particles at any time t , V (t ) = V I (t ) + V N I (t ), can be estimated as follows: where V (t ) denotes the random variable for the number of virus particles at time t , ϕ(t ) the survival probability of the branching process until time t and E[V (t );V (t ) > 0] the expectation of V (t ) for a surviving trajectory until time t . Similar arguments have also been used in the models for HIV (Konrad et al., 2017;Conway and Perelson, 2018). To compute the time for the viral load to reach a certain threshold we set ϕ(t ) = ϕ. In other words, we approximate the survival of the branching process until time t by the total establishment probability expressed in Eq. (4) in the main text. This is a good approximation if the 'typical' time t to reach the threshold is large enough, so that ϕ(t ) is already close to the limit survival probability ϕ. The other term on the right-hand side in Eq. (F.1), the mean of the stochastic process E[V (t )], can be approximated by the deterministic model of the within-host model defined in Eq. (1) in the main text.
As explained in the main text, we set the threshold viral load to 2, 000 virions (Fig 4 in the main text). The mean time to reach this threshold value is then approximated by the time when the size 2, 000 × ϕ is reached in the deterministic model.

F. . Growth rate of the viral population to leading order
The exponential growth rate of the deterministic model described in Eq. (1) in the main text is given by the leading eigenvalue of the system when evaluated at the origin, i.e. at zero virions Bonhoeffer et al. (1997). For efficacies close to the critical efficacy, the eigenvalue is small and can therefore be approximated by the root of a linear equation instead of a higher order polynomial. This approximation yields Under treatment that increases the infected cell death δ, establishing trajectories reach the detectable viral load almost immediately. In contrast, drugs that directly affect the number of virus, i.e. clearance c or production p, allow for trajectories that fluctuate much more, explaining the larger average detection times and the larger variation of detection times for these scenarios.

F. . Explaining the shape of the curves in Fig of the main text
In this section, we provide more detailed explanations about the shapes of the establishment time curves depending on the mode of action of the drug. Throughout this discussion, it is important to keep in mind that for the average establishment times, only trajectories that result in establishment are taken into account. To ease the discussion, Fig D shows Fig 4 from the main text. Treatment that targets the virus infectivity β does not increase the establishment time because these drugs do not affect the virus dynamics itself. Conditioned on virus establishment, the initially present virus particle will infect a target cell relatively quickly, i.e., on a similar time scale than without treatment, and then follow the same dynamics as without treatment. Since the burst size largely exceeds the detection threshold, in our model a single infected cell is sufficient to reach this threshold. Therefore, the establishment time remains largely unaffected by drugs targeting the infectivity β.
For drugs increasing the infected cell death rate δ, the trajectories that contribute to the results in Fig D are the ones that produced a large number of virus particles from a single cell in a short time. This is because of the strongly increased cell death rate for large values of efficacy ε δ . Therefore, a surviving virus trajectory needs to reach large numbers of virus particles in a short time to avoid extinction. This is different for a reduced viral production p where the infected cell death rate is unaffected. Therefore, it is not necessary for a surviving virus trajectory to reach high viral loads very quickly, even though this is of course possible which is reflected by the large 90% confidence interval. This is visualized in Fig D, panel B: green trajectories correspond to drugs affecting the cell death rate and blue trajectories correspond to drugs reducing viral production.
Lastly, increasing the viral clearance rate c by prophylactic treatment increases the establishment time with increasing efficacy, but not as much as treatment with drugs that reduce viral production p. The reason here is that clearance acts just after the viral production, i.e., there is time passing between the production of a virus particle and its clearance. Hence, reducing virus production has a stronger effect on the establishment time than an increase of viral clearance c which acts later in the viral life cycle.

G. Parameter estimation
Patient data from Young et al. (2020) were fitted using the set of differential equations presented in Eq. (1) in the main text. To ensure identifiability of critical parameters of the viral dynamics, i.e. the basic reproductive number R 0 , the loss rate of infected cells δ and the viral production p, the remaining parameters c, k and V 0 were fixed. Viral clearance c was fixed to 10 day −1 . For the eclipse phase k we chose 5 day −1 and the initial inoculum V 0 was set to 1/30 copies.mL −1 (see Gonçalves et al. (2020) for further details). Parameters were estimated in a non-linear mixed effect model using the SAEM algorithm implemented in Monolix (www.lixoft.com). The best fit using all available patient data resulted in the parameter values R 0 = 7.69, δ = 0.595 and p = 11, 200, the principal data set used in the main text (LowN parameter set).

G. . Parameter estimates for individuals plotted in Fig in the main text
Applying this method to data from the 13 untreated patients in Young et al. (2020), we obtain the best parameter set for each individual. The individual parameter sets from four patients (patients 2,4,11,18) were used to plot the deterministic curves in Fig 1A in Table S2. The initial amount of virus particles per mL, V I (0) = 1/30, corresponds to 1 infectious virus particle in absolute numbers in the total upper respiratory tract, which we assume has a volume of 30 mL. The dotted line depicts the detection threshold set to 10 1.84 .

G. . Sensitivity analysis with respect to variations in the fraction of infectious virus particles η
We evaluate how different choices of η, the fraction of infectious virus among all produced virus particles, affect the estimates of the within-host reproductive number R 0 and the burst size η. In the main text, we have used the parameter estimate with η = 10 −3 which resulted in R 0 = 7.69 and N = 18, 823. For a larger fraction of infectious virus particles, η = 10 −2 , we find R 0 = 5.3 and N = 3, 303; for a smaller fraction of infectious virus particles, η = 10 −4 , we obtain R 0 = 9.2 and N = 349, 367. While the within-host reproductive number R 0 does not vary too much between the different choices of η, the burst size N shows large variation. This has no effect on our results on the establishment of a SARS-CoV-2 infection because the burst size always enters in the form of a product with η. In all the different scenarios above, the product η × N varies between 18 for η = 10 −3 and 35 for η = 10 −4 . Overall, the differences in estimates for R 0 will affect the precise estimate of the critical efficacy and differences in the estimate for N translate to differences in the quantitative values of the establishment probability curves below the critical efficacy. The predictions on the detection and extinction time strongly depend on the overall burst size N so that these will vary considerably depending on the choice of η.