Stochastic resetting antiviral therapies prevents drug resistance development
Angelo Marco Ramoso, Juan Antonio Magalang, Daniel Sánchez-Taltavull, Jose Perico Esguerra, Édgar Roldán
SStochastic resetting antiviral therapies prevents drug resistance development
Angelo Marco Ramoso , Juan Antonio Magalang ,Daniel S´anchez-Taltavull , Jose Perico Esguerra , and ´Edgar Rold´an Theoretical Physics Group, National Institute of Physics,University of the Philippines, Diliman, Quezon City 1101, Philippines Department for Visceral Surgery and Medicine,Bern University Hospital, University of Bern, Switzerland ICTP - The Abdus Salam International Centre for Theoretical Physics, Strada Costiera 11, 34151, Trieste, Italy (Dated: September 16, 2020)We study minimal mean-field models of viral drug resistance development in which the efficacyof a therapy is described by a one-dimensional stochastic resetting process with mixed reflecting-absorbing boundary conditions. We derive analytical expressions for the mean survival time for thevirus to develop complete resistance to the drug. We show that the optimal therapy resetting ratesthat achieve a minimum and maximum mean survival times undergo a second and first-order phasetransition-like behaviour as a function of the therapy efficacy drift. We illustrate our results withsimulations of a population-dynamics model of HIV-1 infection.
Antiviral and antiretroviral therapies are continuouslybeing developed to tackle viral diseases [1]. Because ofviral evolution, a therapy that is effective today may notremain effective forever [2]. This process is known as drugresistance development and it is of special importance inchronic infections such as HIV-1 [3–5], and also for herpesand hepatitis [6]. Viral evolution results into an increasein the number of infected cells, leading to a failure of thenormal functions of the infected patient that can resultin their death [2]. Because viruses replicate and mutaterapidly, and mutation is a random process, viral evolu-tion is intrinsically noisy [7, 8]. When possible, prac-titioners overcome this situation by changing the ther-apy given to the patient. Changes in antiviral therapiescan occur because of a detected viral resistance [3, 9],or for purely stochastic reasons e.g. the appearance inthe market of a therapy that is cheaper or has less sec-ondary effects. It remains an open yet challenging prob-lem to characterize and design therapy change protocolsthat ensure large patient survival times that are robustto drug-resistance fluctuations.Stochastic models have been used to study viral evo-lution [7–11]. Examples include the usage of branchingprocess to describe viral persistence and extinction [7, 8],and population dynamics models to assess the impact of atreatment in drug resistance development in the contextof HIV-1 [12]. A model that described RNA virus evolu-tion as a diffusion process in a fitness landscape [13] wasable to explain the experimentally-observed rapid initialgrowth followed by a slower stage of linear growth in thelogarithm of fitness of clone colonies of vesicular stom-atitis virus [14, 15]. In the same vein, the notion of afluctuating therapy efficacy can be used to describe theevolution of therapy protocols.Viral evolution can be modelled as a stochastic dif-fusion process involving incremental changes in efficacy—due to viral mutation— punctuated by sudden changesin efficacy —due to changes in therapy—. Diffusions
FIG. 1. Illustration of the stochastic-resetting model of an-tiviral therapy efficacy. Top: Sketch of the main ingredientsof the model. Bottom: Sample trajectory of the therapy ef-ficacy as a function of time, with A-D illustrating key eventsduring a single realization of the process. A) Initial condi-tion, for a 50% efficient therapy η = 1 /
2. B) Change of thecurrent therapy at a random time (stochastic resetting) bya new one with 50% efficacy. C) Reflecting boundary condi-tion modelling the maximum possible therapy efficacy η = 1.D) Absorbing boundary condition yielding to the dead of thepatient due to a completely inefficient therapy η = 0. with stochastic resetting [16–27] thus provide a promis-ing framework to model therapy evolution. This frame-work has been instrumental to describe biophysical pro-cesses with sudden changes, namely RNA polymerasebacktracking [28, 29], receptor dynamics [30], cell crawl-ing [31], and population dynamics [32–34], see [17] for areview of applications. It remains unclear how in mul- a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p tiscale processes stochastic resetting of a slow variable(e.g. a therapy efficacy) affects a fast variable describingthe dynamics of a time-varying population (e.g. cells andviruses).In this Letter, we introduce a stochastic-resettingmodel for the evolution of the efficacy of an antiviral ther-apy, and study its evolution under different treatmentprotocols. We first describe the therapy efficacy as a one-dimensional resetting biased diffusion model with mixedabsorbing and reflecting boundaries, calculate an exactanalytical expression for the mean first passage time, andtest our result by comparing it with Langevin dynamicssimulations. We then discuss the clinical effects of our vi-ral evolution model by coupling the stochastic resettingbiased-diffusion therapy efficacy to a population dynam-ics model of HIV-1 chronic disease.We describe the efficacy of an antiviral treatment as abounded stochastic process 0 ≤ η ( t ) ≤ η = 0 and η = 1 correspond to a completely ineffective and com-pletely effective therapy, respectively. We assume thattreatments that stop working lead to the death of thepatient, i.e. η = 0 is an absorbing boundary. On theother hand η = 1 is a reflecting boundary set at themaximum 100% therapy efficacy. We model the evolu-tion of therapy efficacy η ( t ) as a biased random walk ordrift-diffusion process in ”efficacy space” with diffusioncoefficient D and drift v . The therapy is often biased inthe direction of lesser efficacy, i.e. the drift is negative v <
0, due to the fact that viruses develop resistanceto the existing therapies that survive in the long run.Furthermore, we complement the biased diffusion modelwith a resetting protocol that switches the therapy effi-cacy to a value η instantaneously at random Poissoniantimes with rate r . Such therapy ”resetting” events canbe due to the introduction of a new dose of drug, thediscovery of more effective variants of the therapy, etc.See Fig. 1 for an illustration of the model and a samplestochastic trajectory of the therapy efficacy.The evolution of the model can be described by theFokker-Planck equation with source terms ∂ t P + ∂ η ( vP − D∂ η P ) = − rP + rδ ( η − η ) , (1)where P ≡ P ( η, t | η ,
0) is the conditional probabilitydensity that the therapy is at η at time t given thatits initial value (at which it is reset at rate r ) was η .The dynamics is complemented by a absorbing boundarycondition at η = 0, lim η → + P ( η, t ) = 0 and a zero-fluxreflecting boundary condition at η = 1, J (1 , t ) = 0 where J ( η, t ) = vP ( η, t ) − D∂ η P ( η, t ) is the probability current.In the following, we derive analytical expressions forthe finite-time survival probability and the mean survivaltime. We denote by survival time the first-passage timeof the therapy efficacy to reach the absorbing boundary η = 0. To this aim, we first make use of a relation be-tween the finite-time survival probability with resetting S r ( η , t ) and the survival probability without resetting S ( η , t ) [20] S r ( η , t ) = e − rt S ( η , t )+ r (cid:90) t d τ e − rτ S ( η , τ ) S r ( η , t − τ ) . (2)From Eq. (2) we show that the Laplace transform ofthe first passage probability with resetting ˜ S r ( η , s ) andwithout resetting ˜ S ( η , s ) are related through the iden-tity (see Supplemental Material)˜ S r ( η , s ) = ( r + s ) ˜ S ( η , s + r ) s + r ˜ S ( η , s + r ) . (3)We show that the Laplace transform of the first passagesurvival probability without resetting is given by˜ S ( η , s ) = e − η v D Dω cosh[ ω ( η − v sinh[ ω ( η − Dω cosh ω − v sinh ω (4)where ω = √ v + 4 Ds/ D . SubstitutingEq. (4) in Eq. (2) and using the relation (cid:104) τ r (cid:105) = − lim s → ∂ ˜ S r ( η , s ) /∂s we obtain the follow-ing analytical expression for the mean first-passage time (cid:104) τ r (cid:105) = φ (Pe , Ω , η ) r (cid:110) Pe sinh (cid:104) Ω (cid:0) η − (cid:1)(cid:105) + Ω cosh (cid:104) Ω (cid:0) η − (cid:1)(cid:105)(cid:111) , (5)where the non-trivial function φ (Pe , Ω , η ) =Pe { e Pe η cosh[Ω( η − − cosh[2Ω( η − } +PeΩ { e Pe η sinh[Ω (cid:0) η − (cid:1) ] − sinh[2Ω( η − } + e Pe η cosh[Ω η ] − cosh[2Ω( η − e Pe η cosh[Ω( η − − v/ D, Ω = √ v + 4 Dr/ D ,and ξ = r/ D .Figure 2 shows an excellent agreement between Eq. (5)and numerical Langevin-dynamics simulations of themodel, for different parameter values. For positive val-ues of v , the mean survival time decreases monotonouslywith the therapy resetting rate, hence it is beneficial toswitch slowly among beneficial therapies (i.e. keeping r small) in order to maximize the survival time. For smalland even negative values of the efficacy drift, the meanfirst-passage time is non monotonous; the minimum aver-age survival time takes place for intermediate values of r .For therapies with large negative bias v <
0, a case thatis is relevant in the context of viral evolution, (cid:104) τ r (cid:105) in-creases monotonically with r , i.e. the maximum averagesurvival is achieved switching the therapies as frequentas possible.In a real-world scenario, practitioners have to deal withlimited resources such as a finite number of therapies dur-ing the life of a patient. Within this scenario, it is im-portant to know what is the optimal resetting rate thatachieves a desired value of the mean survival time of thepatient. To study this problem, we evaluate in Fig. 3 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼ ▼■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ ◆ r < t r > FIG. 2. Mean survival time for a therapy to become com-pletely inefficient, as a function of the therapy resetting rate:numerical simulations (symbols) and analytical results (lines)given by Eq. (5). The data corresponds to parameter valuesdiffusion D = 1 .
5, initial efficacy η = 1 / v : v = 2 ( (cid:7) ), v = 1 ( (cid:4) ) v = 0 ( (cid:72) ), v = − (cid:78) ), and v = − • ). For all parameter values, the sim-ulations were done using Euler’s numerical integration withparameters: 10 number of simulations, each with time step∆ t = 10 − and total simulation time t sim = 10 . the resetting rates r min and r max for which the meansurvival time attains respectively its minimum and max-imum value, within a finite range of resetting rate. Wefind that r min = 0 for rapidly evolving virus ( v large andnegative), whereas r min increases monotonously with v for larger values. Notably r min presents a non-analyicbehaviour at a critical value of v , and hence of the Pecletnumber, as reported by recent work on drift-diffusion pro-cesses with one absorbing boundary [26]. On the otherhand, our results show that the ”optimal” resetting rate r max achieving the maximum mean survival time displaysa dependency with the therapy drift that has reminis-cences of a first-order phase transition; it exhibits a sud-den jump from the maximum allowed resetting rate (forrapidly evolving virus) to zero. Such transition occurs ata critical value of v that depends on the fluctuations D of the therapy efficacy; it can take place at biologicallyrelevant values ( v negative) when D is large enough.We now consider a stochastic mean-field population-dynamics model describing a multicellular organism con-taining healthy H , latent L , and productively infective I cells by a chronic viral disease such as HIV-1. Thestate of the system is described by its time-dependentnumbers H , I and L , whose dynamics is driven by thestochastic-resetting therapy efficacy η . We remark thatwe include in the model a population of latent cells to ac-count for a chronic infection, inspired in previous mathe-matical models of HIV-1 [4]. The dynamics of the modelis given by three coupled ordinary differential equations[Eqs. (6-8) below] driven by an autonomous stochastic ● ● ● ● ● ● ● ● ● ● ●▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ - - v r m i n ● ● ● ● ● ● ● ● ● ● ●▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲ ▲■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ - - v r m ax FIG. 3. Comparison between simulations (symbols) and the-oretical results (lines) of optimal resetting rates as a functionof the therapy efficacy drift: r min (top) and r max (bottom) de-note the value of the resetting rate at which the mean survivaltime attains its minimum and maximum values, respectively.The different curves and symbols are obtained by imposing amaximum allowed resetting rate of 15, for different values of v and D = 0 . • , red line), D = 1 ( (cid:78) , green line) and D = 1 . (cid:4) , purple line). The rest of the simulation parameters areset to the same values as in Fig. 2. differential equation [Eq. (9) below]:d H d t = α − λ H H − (1 − η ) βHI (6)d L d t = (cid:15) (1 − η ) βHI + pL (cid:18) − LK (cid:19) − a L L − λ L L (7)d I d t = (1 − (cid:15) )(1 − η ) βHI + a L L − λ I I (8)d η = (1 − χ )( v d t + √ D d W ) + χ ( η − η ) . (9)Here, α denotes the rate of recruitment of new healthycells, λ H is the death rate of healthy cells, β is the in-fection rate, (cid:15) is the probability of an infection resultinginto a latent cell, p is the proliferation rate of latent in-fected cells, K is the carrying capacity which introducesa logistic growth, a L the activation rate of latent cells, λ L death rate of latent cells, and λ I the death rate of in-fected cells. Note that the infection rate β is multipliedby the instantaneous probability (1 − η ) of the infection tooccur. Finally, χ is a binary variable which equals to one(zero) when a reset occurs (does not occur) with proba-bility r d t (1 − r d t ), and W is the Wiener process. Hence,Eq. (9) complemented with mixed absorbing boundaryconditions describes the previously introduced therapyefficacy stochastic model. Further details of the modelcan be found in the Supplemental Material. M illi on s o f hea l t h y c e ll s η TherapyresetFraction healthy cells η Time (years) 0 0.2 0.4 0.6 0.8 1
Fraction healthy cells η FIG. 4. Numerical simulations of the viral evolutionpopulation-dynamics model given by Eqs. (6-9) in the ab-sence (top) and in the presence (bottom) of therapy re-setting: fraction of healthy cells with respect to the ini-tial value
H/H (green line) and therapy efficacy (blackline) as a function of time. The inset in the top panelshows the average number of healthy cells as a function ofthe therapy efficacy in the absence of resetting, and thered dashed line in the bottom panel illustrates the time atwhich the antiviral therapy efficacy is restored to its initialvalue. Parameters of the simulations: v = − × − days − and D = 10 − days − , α = 6000 days − ml − , λ H =0 .
01 days − , β = 5 × − ml days − , (cid:15) = 0 . , p =0 . − , a L = 0 . − , λ L = 0 .
01 days − , λ I =1 days − , K = 100 cells ml − , with r = 0 (top) and r =(1 / − (bottom), with initial condition η = 0 . H =6 × cells ml − , L = 1 cells ml − , I = 0 cells ml − , andstep size ∆ t = 1days. The inset (up) shows the value of H atthe non-trivial fixed point of the system Eqs. (6-8) for a fixedvalue of η . Next, we illustrate the model with numerical simula-tions showing the impact of the therapy efficacy in thenumber of healthy cells of a patient. For this purpose wenumerically integrate Eqs. (6-9). In this case, we onlyconsider therapies with v < η through the absorbing barrier, followed by the healthycells (Fig. 4, top). Under stochastic changes in the ther-apy η still drifts through the absorbing state barrier, how-ever, after the stochastic reset, we often observe a periodin which the healthy cells recover, delaying the absorp- -0.6 -0.4 -0.2 v (years -1 )
246 -0.6 -0.4 -0.2 v (years -1 ) D ( - y ea r s - ) D ( - y ea r s - ) A BC D Survival time Patient survival time FIG. 5. Mean survival times (in years) as a function of thetherapy efficacy drift v and diffusion D . (A,C) Analyticalvalue of the mean first-passage time for the therapy efficacy η to reach the absorbing boundary η = 0, given by Eq. (5).(B,D) Numerical value of the mean time elapsed until thefraction of healthy cells falls below 1 / r = 0 (A,B) and r = (1 / − (C,D). In (B,D) the valuesof the simulation parameters were set to the same values as inFig. 4 except the time step ∆ t = 0 . tion time of η (Fig. 4, bottom). We remark that, when η reaches the absorbing boundary at zero, the cell pop-ulation still evolves, towards its fixed point. Note that,unlike in our first model, we allow here for resets from η = 0 to η because even in the absence of therapy thepatient can survive until their therapy changes. This mo-tivates us to study first-passage times in the context ofhealthy cells. In doing so, we define the patient survivaltime as the time elapsed until H ≤ α/ (2 λ H ), which cor-responds to the time until it falls below half of the fixedpoint in the absence of infection.We now determine the impact of changes in the ther-apy in a region of parameters for v and D . Figure 5shows the mean survival time (left panels) and the meanpatient survival time (right panels) as a function of v and D obtained from analytical and numerical calcula-tions. For the parameter values studied here, resettingtherapies increases both the mean survival time and themean patient survival time. Their qualitative behaviouris similar: when the drug resistance develops slowly ( v negative but small) the survival time and the patientsurvival time are large, and vice versa. The larger theresistance fluctuations D , the lower the survival times,however this dependency is weaker than for v . Notably,when executing therapy resets at a rate of (1 / − ,the mean patient survival time can exceed 20 years forsmall values of v and D (Fig. 5D).We have introduced a multi-scale stochastic modellinking drug resistance development described by a one-dimensional stochastic resetting process with cell pop-ulation dynamics. We have found that the number ofhealthy cells in a chronic disease such as HIV-1 displaynegative correlation with the drug resistance. Our ana-lytical and numerical results quantify the beneficial as-pects of therapy changes at random Poissonian times forthe mean survival time of a patient as a function of the vi-ral evolution parameters. We have derived an analyticalexpression for the mean survival time of the stochastic-resetting therapy efficacy as a function of its drift and dif-fusivity. This expression can be used to estimate patientsurvival times in some limits. It will be interesting to ex-tend our work to e.g. compare different therapy resettingprotocols under time constraints, account for drug resis-tance development which depends on the instantaneousviral load, account for multiple therapies [35], etc. We ex-pect potential applications of our work to determine theepidemiological impact of drug resistance development,using models where the immunity to certain drugs canbe transported between individuals. [1] M. Lichterfeld and K. C. Zachary, Ther. Adv. ChronicDis. , 293 (2011).[2] D. S. Clutter, M. R. Jordan, S. Bertagnolio, and R. W.Shafer, Infect. Genet. Evol. , 292 (2016).[3] D. Pillay and M. Zambon, Br. Med. J. , 660 (1998).[4] A. S. Perelson and P. W. Nelson, SIAM Rev. , 3 (1999).[5] D. Sanchez-Taltavull, A. Vieiro, and T. Alarcon, J. Math.Biol. , 919 (2016).[6] L. Strasfeld and S. Chou, Infectious Disease Clinics ,809 (2010).[7] L. G. Fabreti, D. Castro, B. Gorzoni, L. M. R. Janini,and F. Antoneli, Bulletin of mathematical biology ,1031 (2019).[8] S. C. Manrubia and E. L´azaro, Phys. Life Rev. , 65(2006).[9] A. M. Wensing, V. Calvez, F. Ceccherini-Silberstein,C. Charpentier, H. F. G¨unthard, R. Paredes, R. W.Shafer, and D. D. Richman, Top. Antivir. Med. , 111(2019).[10] F. Tria, M. Laessig, L. Peliti, and S. Franz, J. Stat. Mech. , P07008 (2005).[11] M. I. Nelson, L. Simonsen, C. Viboud, M. A. Miller,J. Taylor, K. St George, S. B. Griesemer, E. Ghedin,N. A. Sengamalay, D. J. Spiro, et al., PLoS Pathog. ,e125 (2006).[12] C. Zitzmann and L. Kaderali, Front. Microbiol. , 1546(2018).[13] L. S. Tsimring, H. Levine, and D. A. Kessler, Phys. Rev.Lett. , 4440 (1996).[14] I. S. Novella, E. A. Duarte, S. F. Elena, A. Moya,E. Domingo, and J. J. Holland, PNAS , 5841 (1995). [15] J. J. Holland, J. C. De La Torre, D. Clarke, andE. Duarte, J. Virol. , 2960 (1991).[16] M. R. Evans and S. N. Majumdar, Phys. Rev. Lett. ,160601 (2011).[17] M. R. Evans, S. N. Majumdar, and G. Schehr, J. Phys.A , 193001 (2020).[18] U. Bhat, C. De Bacco, and S. Redner, J. Stat. Mech. , 083401 (2016).[19] (cid:32)L. Ku´smierz and E. Gudowska-Nowak, Phys. Rev. E ,052127 (2015).[20] A. Pal and V. Prasad, Phys. Rev. E , 032123 (2019).[21] A. Pal and V. Prasad, Phys. Rev. Res. , 032001(R)(2019).[22] X. Durang, S. Lee, L. Lizana, and J.-y. Jeon, J. Phys. A , 224001 (2019).[23] C. Christou and A. Schadschneider, J. Phys. A ,285003 (2015).[24] A. Pal and S. Reuveni, Phys. Rev. Lett. , 030603(2017).[25] A. Chechkin and I. Sokolov, Phys. Rev. Lett. , 042128(2019).[26] S. Ray, D. Mondal, and S. Reuveni, J. Phys. A , 255002(2019).[27] U. Basu, A. Kundu, and A. Pal, Phys. Rev. E ,032136 (2019).[28] ´E. Rold´an, A. Lisica, D. S´anchez-Taltavull, and S. W.Grill, Phys. Rev. E , 062411 (2016).[29] G. Tucci, A. Gambassi, S. Gupta, and ´E. Rold´an, arXivpreprint arXiv:2005.05173 (2020).[30] T. Mora, Phys. Rev. Lett. , 038102 (2015).[31] P. C. Bressloff, Phys. Rev. E , 022134 (2020).[32] G. Mercado-V´asquez and D. Boyer, J. Phys. A ,405601 (2018).[33] T. T. da Silva and M. D. Fragoso, J. Phys. A , 505002(2018).[34] R. Garc´ıa-Garc´ıa, A. Genthon, and D. Lacoste, Phys.Rev. E , 042413 (2019).[35] D. D. Richman, D. M. Margolis, M. Delaney, W. C.Greene, D. Hazuda, and R. J. Pomerantz, Science ,1304 (2009).[36] S. Redner, A guide to first-passage processes (CambridgeUniversity Press, 2001).[37] ´E. Rold´an and S. Gupta, Phys. Rev. E , 022130 (2017).[38] T. Pierson, J. McArthur, and R. F. Siliciano, Annu. Rev.Immunol. , 665 (2000). SUPPLEMENTAL MATERIALS1. Analytical expression for the mean survival time
In this section we provide additional details about the derivation of Eq. (5) in the Main Text for the mean absorptiontime for the one-dimensional stochastic-resetting process that we use to describe the therapy efficacy. The derivationproceeds as follows: First, we derive an analytical expression for the survival probability in the absence of resetting r = 0. Next, we use a known relation to derive the Laplace transform of survival probability with resetting from thesurvival probability without resetting. Third, we derive the mean survival time from the Laplace transform of thesurvival probability with resetting. First passage without resetting
We first consider for the ease of analytical calculations the case r = 0 in which the therapy efficacy is a drift diffusionprocess. We denote by P ( η, t ) the conditional probability density that the efficacy of the therapy is η at time t , giventhat its initial value at time t = 0 was η . It obeys the Fokker-Planck equation that results from taking r = 0 inEq. (1) in the Main Text ∂P ( η, t ) ∂t + ∂J ( η, t ) ∂η = 0 , (S10)where J ( η, t ) is the probability current defined in this case as J ( η, t ) ≡ vP ( η, t ) − D ∂P ( η, t ) ∂η . (S11)Equations (S10) and (S11) are complemented with the initial and boundary conditions of the process P ( η, t = 0) = δ ( η − η ) , (S12)lim η → + P ( η, t ) = 0 , (S13) J ( η = 1 , t ) = 0 . (S14)Here, Eq. (S12) accounts for the initial condition η (0) = η , Eq. (S13) for the absorbing boundary at η = 0, andEq. (S14) for the reflecting boundary condition at η = 1. Taking the Laplace transform on Eq. (S10), we get s P ( η, s ) − P ( η, t = 0) + v P (cid:48) ( η, s ) − D P (cid:48)(cid:48) ( η, s ) = 0 , (S15)where the prime denotes derivative with respect to η and we have introduced the notation P ( η, s ) ≡ (cid:90) ∞ e − st P ( η, t ) d t, (S16)for the Laplace transform of the propagator. Applying the initial condition (S12) into Eq. (S15), we obtain s P ( η, s ) − δ ( η − η ) + v P (cid:48) ( η, s ) − D P (cid:48)(cid:48) ( η, s ) = 0 . (S17)We now concentrate our efforts in obtaining an analytical solution to Eq. (S17). First we express P ( η, s ) as apiecewise continuous function: P ( η, s ) = (cid:40) P < ( η, s ) when η < η P > ( η, s ) when η > η , (S18)for η (cid:54) = η . Both P < and P > obey Eq. (S17) in η < η and η > η respectively, i.e. D P (cid:48)(cid:48) < ( η, s ) − v P (cid:48) < ( η, s ) − s P < ( η, s ) = 0 , (S19) D P (cid:48)(cid:48) < ( η, s ) − v P (cid:48) < ( η, s ) − s P < ( η, s ) = 0 . (S20)The solutions of Eqs. (S19-S20) are given by sums of exponential functions P < ( η, s ) = a + e α + η + a − e α − η , (S21) P > ( η, s ) = b + e α + η + b − e α − η . (S22)Plugging in Eq. (S21) in (S19) and Eq. (S22) in (S20) yields the second order equation Dα ± − vα ± + s = 0 , (S23)whose solutions are given by α ± = v D ± w, (S24)where w ≡ √ v + 4 Ds D . (S25)The parameters a + , a − , b + , and b − can be determined using the boundary conditions. Applying the absorbingboundary condition (S13) into Eq. (S21) we get a + + a − = 0, hence defining a ≡ a + , we find that P < has the form P < ( η, s ) = a (cid:104) e α + η − e α − η (cid:105) . (S26)Next, we use the reflecting boundary condition (S14) in (S22) which implies v P > ( η = 1 , s ) = D P (cid:48) > ( η = 1 , s ). Thisresults in the following relation v (cid:104) b + e α + + b − e α − (cid:105) = D (cid:104) α + b + e α + + α − b − e α − (cid:105) , or equivalently ( v − Dα + ) e α + b + =( Dα − − v ) e α − b − ≡ b . Thus we get P > ( η, s ) = b (cid:20) e α + ( η − v − Dα + − e α − ( η − v − Dα − (cid:21) . (S27)The continuity condition of P ( η, s ) at η = η in Eqs. (S26-S27) implies that b (cid:20) e α + ( η − v − Dα + − e α − ( η − v − Dα − (cid:21) = a [ e α + η − e α − η ] ≡ c (S28)therefore P < ( η, s ) = c (cid:20) e α + η − e α − η e α + η − e α − η (cid:21) , (S29) P > ( η, s ) = c e α + ( η − v − Dα + − e α − ( η − v − Dα − e α + ( η − v − Dα + − e α − ( η − v − Dα − . (S30)To solve for c , we solve the differential equation (S17) at η = η , D P (cid:48)(cid:48) ( η, s ) = v P (cid:48) ( η, s ) + s P ( η, s ) − δ ( η − η ) . (S31)Take the integral both sides with respect to η from η − (cid:15) to η + (cid:15) and then the limit of (cid:15) to 0 we obtain D (cid:104) P (cid:48) > ( η , s ) − P (cid:48) < ( η , s ) (cid:105) = v (cid:104) P > ( η , s ) − P < ( η , s ) (cid:105) − P (cid:48) > ( η , s ) − P (cid:48) < ( η , s ) = − D . (S33)Using Eqs. (S29) and (S30) in (S33) and solving for c we get, after some cumbersome simplifications, c = 2 (cid:32) √ v + 4 Ds cosh (cid:34) √ v + 4 Ds ( η − D (cid:35) + v sinh (cid:34) √ v + 4 Ds ( η − D (cid:35)(cid:33) sinh (cid:34) √ v + 4 Dsη D (cid:35) ( v + 4 Ds ) cosh (cid:34) √ v + 4 Ds D (cid:35) − v √ v + 4 Ds sinh (cid:34) √ v + 4 Ds D (cid:35) . (S34)The Laplace transform of the first-passage (survival) time probability at η = 0 without resetting is given by [36] F ( s ) = D ∂∂η P < ( η, s ) (cid:12)(cid:12)(cid:12)(cid:12) η =0 − v P < ( η, s ) (cid:12)(cid:12)(cid:12)(cid:12) η =0 , (S35)where the second term at the right-hand side vanishes because of the absorbing boundary condition. Note that herewe have introduced the subscript 0 to emphasize that this first-passage statistic refers to the case r = 0. Since c doesnot depend on η , Eq. (S35) together with Eq. (S29) imply that F ( s ) = cD α + − α − e α + η − e α − η . (S36)Substituing the values of α + and α − [Eq. (S24)] and C [Eq. (S34)] in Eq. (S36) we arrive at the analytical expressionfor the Laplace transform of the first-passage-time probability for the drift-diffusion process with mixed boundaryconditions: F ( s ) = √ v + 4 Ds (cid:32) √ v + 4 Ds cosh (cid:34) √ v + 4 Ds ( η − D (cid:35) + v sinh (cid:34) √ v + 4 Ds ( η − D (cid:35)(cid:33) exp (cid:104) η v D (cid:105) (cid:32) ( v + 4 Ds ) cosh (cid:34) √ v + 4 Ds D (cid:35) − v √ v + 4 Ds sinh (cid:34) √ v + 4 Ds D (cid:35)(cid:33) . (S37)From Eq. (S37) we can also derive an analytical expression for the Laplace transform of the survival probability,which is defined as follows S ( s ) ≡ (cid:90) + d η P ( η, s ) (S38)= (cid:90) η + P < ( η, s ) d η + (cid:90) η P > ( η, s ) d η, (S39)where in the second line we have used Eq. (S18). We now recall the relation (valid for any η >
0) between the survivalprobability and the first passage probability F ( t ) = − ∂ t S ( t ) which takes the form F ( s ) = − [ s S ( s ) − S ( t = 0)]upon taking the Laplace transform, or equivalently F ( s ) = 1 − s S ( s ) . (S40)Here we haved used the fact that at time t = 0, the process has not been absorbed by the absorbing boundary η = 0,i.e. η > First passage with resetting
Following [20, 37], the survival probability with resetting (denoted here with subscript r ) obeys a renewal equationtogether with the survival probability without resetting S r ( t ) = e − rt S ( t ) + r (cid:90) t d τ e − rτ S ( τ ) S r ( t − τ ) (S41)Taking the Laplace transform on both sides of Eq. (S41), we get S r ( s ) = S ( s + r ) + r S ( s + r ) S r ( s ), which solvingfor S r ( s ) yields S r ( s ) = S ( s + r )1 − r S ( s + r ) . (S42)Therefore, the Laplace transform of the first-passage-time probability with resetting can be expressed in terms of theLaplace transform of the survival probability without resetting as follows F r ( s ) = 1 − s S ( s + r )1 − r S ( s + r ) . (S43)Shifting the relation (S40) from s to s + r , i.e. using S ( s + r ) = (1 − F ( s + r )) / ( s + r ) in Eq. (S43) we obtain aftersome simplifications F r ( s ) = (cid:16) r + s (cid:17) F ( s + r ) s + r F ( s + r ) (S44)After some algebra, in using Eq. (S37) in (S44) and the relation (cid:104) τ r (cid:105) = − lim s → ∂ s F r ( s ) we derive Eq. (5) in theMain Text for the mean first-passage time, copied here for convenience (cid:10) τ r (cid:11) = Pe exp (cid:104) Pe η (cid:105) cosh (cid:104) Ω (cid:0) η − (cid:1)(cid:105) − cosh (cid:104) (cid:0) η − (cid:1)(cid:105) +Pe Ω exp (cid:104) Pe η (cid:105) sinh (cid:104) Ω (cid:0) η − (cid:1)(cid:105) − sinh (cid:104) (cid:0) η − (cid:1)(cid:105) − ξ
1+ cosh (cid:104) (cid:0) η − (cid:1)(cid:105) − exp (cid:104) Pe η (cid:105) cosh (cid:104) Ω η (cid:105) − exp (cid:104) Pe η (cid:105) cosh (cid:104) Ω (cid:0) η − (cid:1)(cid:105) r (cid:110) Pe sinh (cid:104) Ω (cid:0) η − (cid:1)(cid:105) + Ω cosh (cid:104) Ω (cid:0) η − (cid:1)(cid:105)(cid:111) (S45)where we have introduced the variablesΩ = √ v + 4 Dr D ; Pe = v D ; ξ = r D . (S46)
S2. Details of the biophysical model
Our biological model [Eqs. (6-9) in the Main Text] is an adaptation of that reported in Ref. [4] and describes theevolution of three different cell populations in blood: healthy cells H which are susceptible to be infected, infectedcells L which are in a latent state, and productively infected cells I which can infect healthy cells. The dynamics isdescribed by the reactions described below, where the chronic feature is modelled by adding a latent population witha logistic growth. Such dynamics have been hypothesized in [38] and used to model HIV-1 dynamics in [4, 5]. • Recruitment of a new healthy cell by the organism at a rate a : ∅ → H . • Death of a healthy at a rate λ H : H → ∅ . • A healthy cell being infected, which we assume it is proportional to the number of infected cells, I , and inverselyproportional to the efficacy of the therapy η , at a rate β . With probability (cid:15) the resulting cell will be a latentcell: H → L , and with probability 1 − (cid:15) the resulting cell will be a productively infected cell: H → I . • Proliferation of a latent cell with rate p : L → L + L . • Death of latent cells by competition for resources, with rate pK : L + L → ∅ . This is an artificial reaction tomimic the logistic growth, see [5] for details in HIV-1 modeling. • Activation of a latent cell resulting into a productively infected cell at rate a L : L → I . • Death of a latent cell at rate λ L : L → ∅ . • Death of a productively infected cell at rate λ I : I → ∅→ ∅