Front propagation speeds of T7 virus mutants

Front propagation speeds of T7 virus mutants V. L. de Rioja,1, J. Fort,1and N. Isern1 1Complex Systems Laboratory, Departament de Física, Universitat de Girona, 17071 Girona, Catalonia, Spain. We propose a new reaction-di¤usion model with an eclipse time to study the spread of viruses on bacterial populations. This new model is both biologically and physically sound, unlike previous ones. We determine important parameter values from experimental data, such as the one-step growth. We verify the proposed model by comparing theoretical and experimental data of the front propagation speed for several T7 virus strains. I. INTRODUCTION Bacterial viruses or bacteriophages (literally ‘eaters of bacteria’) infect and replicate within bacteria. Right after their discovery, phages were used as an early form of biotechnology to …ght bacterial pathogens. Nowadays, drug-resistant strains for many bacteria have appeared and this has lead to a revived interest in this kind of therapy [1]. Moreover, these viruses are among the most common and diverse entities in the biosphere, so it is important to attain a better and more accurate knowledge of their dynamics. Understanding the speed of virus infection fronts is also important in the context of cancer treatment [2]. It is possible to see with the naked eye how the spreading dynamics of viruses works in a medium of susceptible host bacteria. When a small quantity of phages is inoculated into a tiny, central region of liquid agar with host cells (bacteria in our case), the continuous replication and di¤usion of viruses lead to an enlarging dark region, composed of dead cells. Such a region of lysed (i.e. dead) cells, surrounded by unlysed cells, is called a plaque. The growth process starts when a free virus di¤uses into a host bacterium, adsorbs on its surface, injects its DNA into it, replicates within and …nally (after a certain time) the bacterium dies and expels a new generation of viruses. The progeny viruses di¤use to surrounding host cells, and the cycle repeats again. The propagating front has a well-characterized speed, typically less than a millimeter per hour, which has been measured experimentally, and for which we try below to get a realistic and accurate reaction-di¤usion model. Numerous models of phage plaque enlargement have been proposed. The oldest and simplest one is due to Koch [3], who suggested that the di¤usion speed was proportional to qD , where Dis the di¤usion coe¢ cient and the phage latent period (i.e., the time during which bacteriophages are inside cells, and thus not moving). By incorporating additional kinetic parameters, Yin and McCaskill [4] constructed a reaction-di¤usion system and obtained the speed of travelling-wave solutions. Later Electronic address: [email protected] You and Yin [5] supported the previous idea of an existing travelling-wave solution through numerical simulations of the same problem. However the models due to Yin and co-workers [4, 5] lead, for parameter values derived from independent experiments, to speeds much faster than the experimental ones [4, 6]. It was then realized that the delay time or latent period (i.e., the time interval during which a virus is inside a cell and thus does not move) delays virus di¤usion, and that this important e¤ect could explain the slowness of the experimental speeds [7]. By solving the problem numerically, good agreement with experiment was attained (without …tting any parameter values) [7]. However the equations were not fully understood from a biological viewpoint, as we shall explain below. Later Ortega-Cejas et al. [8] obtained some approximate but explicit formulas for the front speed based on the model in Ref. [7]. Amongst more recent models, Amor and Fort [9] proposed a new improved set of equations which satisfactorily explained the observed speeds of VSV (vesicular stomatitis virus) infections, but still with some terms lacking a clear biological interpretation. Using various bacteriophage T7 mutants in a growing plaque on E. coli host bacteria, Yin measured experimentally [10] the radial propagation speed for plaques of three mutant T7 virus strains (namely, p001, p005 and the wild type), …nding di¤erent speeds depending on the type of mutant. These are the experiments that we want to explain. On the basis of the model for VSV infections [9], we rewrite the equations carefully, so that they acquire full biological and mathematical meaning, and we apply them to T7 strains. With this new model we obtain a good agreement with the experimental results in Ref. [10], without requiring the use of any free or adjustable parameters. In this paper, we introduce a new reaction-di¤usion set of equations to explain the existing experimental data on the growth of T7 plaques on bacteria. In section II we present the new time-delayed model and we discuss why our modi…cations are reasonable. Section III is devoted to estimations of the necessary parameter values from independent experiments. In section IV, the results are compared with experimental data for the propagation speed of three strains of the T7 virus, and section V presents a simpli…cation of the model yielding similar results. In Sec. VI we compare to other time-delayed 2 models. Finally, section VII is devoted to …nal conclusions, with particular attention to the model features and how the results are improved over previous models. II. REACTION-DIFFUSION MODEL We model the spatial dynamics of T7 mutants infecting host cells by considering interactions between three species: viruses (V), uninfected bacteria (B) and infected bacteria (I). Those processes can be described schematically as V+Bk1 ! Ik2 ! YV; (1) where k1is the adsorption rate, k2the death rate of infected bacteria, and Y(yield or burst size) is the number of new viruses released per lysed host bacteria. These three parameters (k1,k2and Y) depend on the mutant strain considered. The experiment on which this theoretical work is focused [10] was conducted in agar (so that host bacteria are immobilized) and cells were initially in the stationary phase, i.e. with bacterial growth and death in balance (so that the number density of live bacteria does not change appreciably before viruses arrive). Viruses can move and adsorb on host bacteria, infecting cells and producing new viruses. Some previous models have the drawback of assuming logistic dynamics, namely [7–9] @[I](r; t) @t =k2[I] (r; t)1[I] (r; t) [I]max ;(2) in the absence of uninfected cells (i.e., if all cells are initially infected). [I]in Eq. (2) is the concentration of infected hosts, [I]max their maximum concentration, r the distance from the inoculation point, and tthe time. Let us de…ne ‘free space’ as the fraction of space not occupied by infected cells, relative to the maximum possible value that can be occupied by them, i.e. 1[I] (r; t)=[I]max. Equation (2) describes well the onestep growth experiment (see Fig. 1 in [7]) but has no biological meaning. Indeed, it assumes that the death rate of infected cell is proportional not only to the concentration of infected cells [I](which is reasonable), but also to the free space (term within brackets). Thus, we propose to replace this equation by taking into account the eclipse time between adsorption and the onset of the release of the virus progeny. Therefore, in the absence of adsorption we propose to replace Eq. (2) by @[I](r; t) @t =k2[I] (r; t ):(3) Note that we do not assume that all cells die at the same time after infection. That assumption is made in the perfect delay model [11], which makes use of [V] (r; t ) [B] (r; t )instead of k2[I] (r; t )(as we will see in Sec. VI in detail). But the perfect delay model disagrees with biological experiments (because one-step experiments do not display a vertical step, see Fig. 1). In contrast, Eq. (3) does not represent a perfect vertical step, but a gradual increase after an eclipse time . The model we present below is an alternative to the exponential, non-delayed model [i.e., a term k2[I] (r; t)], and the perfect delay model [i.e., a term [V] (r; t ) [B] (r; t )], and is more realistic than both extreme models. In the presence of adsorption, the model we propose is thus described by (see Appendix A): @[I](r; t) @t =k1[V](r; t)[B](r; t)k2[I](r; t );(4) @[V](r; t) @t + 2 @2[V](r; t) @t2=Deff @2[V](r; t) @r2 +F(r; t) 2k1[V](r; t)@[B](r; t) @t  2k1[B](r; t)F(r; t) +  2k2Y@ @t [[I](r; t )] ;(5) @[B](r; t) @t =k1[V](r; t)[B](r; t);(6) where [V]and [B]are the concentration of viruses and uninfected bacteria respectively, and Deff is the e¤ective di¤usion coe¢ cient of viruses (see next section). Bacteria do not di¤use because they are immobilized by the agar in this experiment. The virus growth function, F(r; t), in Eq. (5) is F(r; t) = k1[V](r; t)[B](r; t) + k2Y[I](r; t ):(7) In this model [Eqs. (4)-(7)], the time derivative @ @t represents the change of the population number over time and the second space derivative @2 @r2is related to the di¤usion through space. Terms proportional to k1account for the decay of viruses [Eqs. (5) and (7)] and host bacteria [Eq. (6)] and the creation of infected cells [Eq. (4)], as a result of the infection process (note that these terms are the same as Eqs. (9) and (11)-(13) in Ref. [9]). Infected cells also decay following their own rate of death k2[Eq. (4)], and as shown by Eq. (1), for each dead cell the viruses increase their number Ytimes [Eqs. (5) and (7)]. The terms proportional to in Eq. (5) are second-order corrections (see Appendix A), they were applied already in Ref. [9], and they take care of the time delay due to the fact that viruses spend a time inside cells before the new generation disperses away [7]. As mentioned above, the main drawback of Ref. [9] is studying the death of the infected cells from a logistic equation, which has no biological sense. Therefore, here we present a new model with two main e¤ects: (i) the second-order correction that has been shown to be fundamental to describe time-delayed biological fronts [7, 9, 12] (i.e., the terms proportional to 3 FIG. 1: One-step growth curves of T7 mutants adapted to the model in this paper. Experimental data (for the wild T7, for the p001 mutant and 4for the p005 mutant) have been obtained from Fig. 3 in Ref. [10]. Full, dashed and dotted lines correspond to the …ts for the wild type and p001 and p005 mutants, respectively. in Eq. (5)), and also (ii) a biologically meaningful description of the death process, Eq. (3) (instead of logistic growth dynamics, Eq. (2)). Other authors have also described the death process through including an eclipse time with terms proportional to concentrations at t, rather than a logistic function [11, 13]. However, those models do not include any second-order terms, i.e. any di¤usive delay (e¤ect (i) in the previous paragraph), which is necessary to take proper account of the fact that viruses do not move during a time interval , because they are inside the infected cells. The death of infected cells is also described in Refs. [11, 13] di¤erently than in our model (in Sec. VI we discuss this in more detail and compare the models and experiments). We introduce dimensionless variables to simplify the analysis. Let B0be the initial concentration of bacteria, then B[B]=B0,V[V]=B0,I[I]=B0,tk2tand rrpk2=Deff are the new dimensionless variables, and k2and k1B0=k2the new dimensionless parameters. The aim is to …nd the speed of travelling-wave solutions which satis…es the set of di¤erential equations (4)-(6). These become single-variable di¤erential equations by using the co-moving coordinate zrct._ c (positive) is the dimensionless wave front speed and is related to the dimensional speed cby _ c=c=pk2Deff . Following previous work [4, 7, 9], we assume that the concentrations at the leading edge of the propagation front (z! 1) are (V ; B; I) = (V;1B; I)'(0;1;0), where ! = (V; B; I) = ! 0exp(z). For non-trivial solutions to exist, the determinant of the matrix corresponding to the linearized model must be zero. This leads us to the following characteristic equation 1 2c2c3+ec c21 +  2ec 2 + 21 +  2Yec ec c +ec  21Y 21= 0:(8) It is known that, according to marginal stability analysis [14], the propagation front moves with the minimum possible speed. Therefore, c= min >0[c()] ;(9) where c()is given implicitly by Eq. (8). III. PARAMETER VALUES We have a new time-delayed model which depends on various parameters. It is necessary to estimate their values from experiments di¤erent from the front-speed experiments that we want to explain. The front propagation speed depends on the viral di¤usivity Deff , the average yield Y, the kinetic parameters k1and k2, the host concentration B0and the eclipse time . Since we aim to explain the experimental data in Ref. [10], these parameters must be determined for strains of the T7 virus and E. coli bacteria. Yin and co-workers noted that the di¤usion coe¢ cient Dof viruses in agar must be corrected by the fact that host bacteria adsorb the viruses, and this leads to more tortuous paths for the viruses at high bacterial concentrations. As noted in previous work, the e¤ective coe¢ cient Deff is therefore given by Fricke’s law [7], Deff =1f 1 + f x D; (10) where fis the initial concentration of bacteria relative to its maximum, i.e. f=B0=Bmax, and xstands as an approximation of the cells’shape. For spherical particles x= 2, while for E. coli it is more accurate to use x= 1:67 [7]. Note that the di¤usivity coe¢ cient Dcorresponds to viruses moving through agar in absence of bacteria (f= 0). T7 viruses are very similar to phage P22 in shape and size, thus we use the corresponding value D= 4 108 cm2/s [7]. The rate of adsorption of viruses, k1, was estimated from a separate experiment conducted in KCN, a substance that prevents viruses from reproducing. We have only one experimental value for the T7 virus, k1= (1:29 0:59)109ml/min [7], corresponding to the wild strains. For the other mutants, we have not been able to …nd any reliable experimental value, thus we will use the 4 same value of k1for all three strains. We will return to discuss this parameter in the next section. Finally, the parameters ,Yand k2are obtained from the so-called one-step growth experiments. They consist in measuring the concentration of viruses as a function of time for a given initial, homogeneous population of infected bacteria. Depending on the T7 mutant, the curves are di¤erent (see Fig. 1) and so will be the parameter values. Fig. 1 allows us to obtain the necessary information from each mutant to estimate its value of ,Yand k2, as we next explain. The eclipse phase of the one-step growth (Fig. 1) corresponds to the stage between adsorption (t= 0) and the …rst release of viruses (i.e., the beginning of the rise in virus density). This interval of time (the eclipse time) is 16 minutes for the wild type and 14 minutes for the p001 and p005 mutants. Note that if we used higher values for ,e.g. = 18 min for the wild strain, for t= 17 min we would have @[I](r;t) @t = 0 according to Eq. (3), whereas we must have @[I](r;t) @t 6= 0 according to Fig. 1. The average burst sizes Ydi¤er signi…cantly for the three mutants. They can be calculated as the quotient between the maximum and the initial concentration of viruses, i.e. Y=Vmax V0, where according to Fig. 1 V0= 2 108ml1for the three kinds of mutants. Inserting the data in Fig. 1, we obtain the yields Y= 34:5 for the wild type, Y= 56:5for the p001mutant, and Y= 65 for the p005 mutant. As we shall see, these higher productivities of new generations for the two mutants result in faster infections (relative to the wild type). The three yields above have been obtained for cells in agar-immobilized microcolonies containing many cells. As noted by Yin and McCaskill [4], such yields are substantially lower than the typical yield for an isolated cell under optimal conditions (Y200). Yin and McCaskill suggested that this di¤erence may be due to a number of factors, such as inherently lower yields per cell when immobilized in agar, premature lysis or inhibition due to the death of adjacent cells, high multiplicities of adsorption required for host infection, readsorption of newly released viruses on cell fragments, etc. [4]. If we used the yield for an isolated cell, we would have to incorporate additional terms to include other possible kinds of death and interactions in our mathematical model. However, the measured experimental values of the burst size quoted above (for cells in agar-immobilized microcolonies containing many cells) implicitly include these possible interactions. The rate of death of infected bacteria k2may be understood as the reproduction of viruses, because viruses replicate as bacteria die. For t <  there are no new viruses (see Fig. 1), so no infected cells have died yet and thus [I] (t) = I0. For tEq. (3) yields d[I] = k2I0dt. Because each infected cell produces Y viruses, d[V] = Y d[I] = k2Y I0dt =k2Vmaxdt. Therefore, the slope of each straight line in Fig. 1 is k2Vmax, and k2=VmaxV0 tVmax 1 t, where tis the time interval during which [V]increases, also known as the rise period 1 k2. It is straightforward to estimate the values of k2from the …gure, and they turn out to be 1=4min1for the wild type and 1=6min1for the p001 and p005 mutants. It is important to remember that all of these parameters are known a priori, thus we do not use any free or adjustable parameters in our predictions. Accordingly to Fig. 1, the average latent period is +1 2k2, where is the eclipse time. For clarity, we mention that when all infected cells have died, no more viruses are produced (Fig. 1, right side) and Eq. (3) obviously breaks down. Thus the general evolution equation we propose is @[I](r; t) @t =k2[I] (r; t )if [V]< Vmax 0if [V] = Vmax (11) where the second line is analogous to some approaches to single-species systems (see Eq. (9) in Ref. [19]). However this point is, in fact, unnecessary for the purposes of the present paper because the front speed is computed at the leading edge of the infection front, where [V]'0(Sec. II). Obviously, in Eq. (11) the condition [V]< Vmax is equivalent to [I]6= 0, and the condition [V] = Vmax is equivalent to [I] = 0. IV. THEORY VERSUS EXPERIMENT In this section we study the spatial dynamics of different T7 virus strains. The experimental data (black squares in Fig. 2) and their error bars were obtained in Ref. [10] for plaques where the concentration of nutrient was 10 g/l, which corresponds to f= 0:2(see Ref. [4], pp. 1543-1544) and Bmax = 107ml1(see Ref. [4], Fig. 3a) thus B0= 2 106ml1. The theoretical results will be calculated below with the parameters Y,k2and for each strain extracted from Fig. 1 (as detailed in Sec. III), and the mean values of k1 and Deff . Because the value of k1is substantially more uncertain than those of other parameters [7], the corresponding error bars are obtained from the experimental range of k1;namely k1= (1:29 0:59) 109ml/min [7]. The classical approach with no delay or eclipse time, due to Yin and McCaskill (triangles in Fig. 2), predicts speeds much faster than the experimental ones (black squares). This model by Yin and McCaskill [4] (with k1= 0, as noted in Ref. [15]) is the same as our model [Eqs. (4)-(7)] with = 0,i.e. @[I](r; t) @t =k1[V](r; t)[B](r; t)k2[I](r; t);(12) @[V](r; t) @t =Deff @2[V](r; t) @r2k1[V](r; t)[B](r; t) +k2[I](r; t)(13) 5 FIG. 2: Front propagation speeds for T7 mutants (wild, p001 and p005). Black squares refer to experimental data and white symbols to the theoretical models: triangles for the classical Yin et al. model, circles for the new model, stars for the simpli…ed model explained in section V, rhombuses for the model by Jones et al. (2012), and white squares for the model by Gourley et al. (2005) both from section VI. @[B](r; t) @t =k1[V](r; t)[B](r; t);(14) The new model introduced in this paper (circles in Fig. 2) agrees better to the experimental data than the classical model Yin and McCaskill, for all three mutants. This improvement is clearly visible in Fig. 2, where we see that the results from the new model lie much closer to the experimental data than the classical model [4]. If we calculate the errors of the models versus the experimental data, the classical model by Yin and McCaskill has an average error of 75%, compared to only 10% for the new model presented here. V. SIMPLIFIED MATHEMATICAL MODEL Our new model, Eqs. (4)-(7), yields a rather complex characteristic equation, Eq. (8), from which we compute the front speeds. In this section we derive a simpli…ed expression leading to similar results. We proceed by removing each term and evaluating its contribution to the front speed, in order to ultimately keep only those terms that have a major contribution on the model results. In this way, it is easy to see that all of the terms in Eqs. (4) and (6) are important to achieve a good result, but some terms in Eq. (5) are not. Hence, we just modify this equation. On one hand, the expansion of F(r; t)to second-order (the three last terms in Eq. (5)) introduces a small change on the results. We can neglect all reaction terms proportional to in this equation. On the other hand, if we understand the right side of Eq. (5) as the di¤usion term, plus the reaction term (plus second-order approximations), we can also neglect the adsorption of virus into bacteria, i.e. the term with k1in Eq. (7). Di¤usion and creation of new viruses are thus the terms with major contributions to the front speed. In this simpli…ed model we can therefore replace Eq. (5) in our set by @[V](r; t) @t + 2 @2[V](r; t) @t2 =Deff @2[V](r; t) @r2+k2Y[I](r; t ):(15) Considering now the set composed by Eqs. (4),(6) and (15) we obtain a new characteristic equation, c +2 2c21c + ec Y ec = 0;(16) much simpler than the previous Eq. (8). The results of this model are shown as stars in Fig. 2. As it can be seen, the front speeds of the simpli…ed model (stars) are always slightly slower than those found by the main model (circles). But the di¤erence between the two models is only about 4% in all three cases. By comparing with the experimental data (black squares in Fig. 2), we see that the simpli…ed model in this section (stars, Eq. (16)) is still much better than the classical one (triangles) in spite of being much simpler than the complete model in Sec. II (circles, Eq. (8)). VI. COMPARISON TO OTHER TIME-DELAYED MODELS Some other authors have also described the death process by considering concentrations at t, rather than a logistic function [11, 13]. However, as mentioned above, those models do not include the di¤usive delay (i.e., second-order corrections), which is necessary because viruses do not di¤use when they are inside the infected cells. Another di¤erence between our model and that in Ref. [11] is that the term k2I(r; t )in our model is replaced by k1B(r; t )V(r; t ). From a conceptual point of view, in our model the infected cells present at the system at time tbegin to die at time t, an do so gradually thereafter (with rate k2). Thus not all cells die exactly at time tin our model, in agreement with the experimental data (Fig. 1). In contrast, according to the model in Ref. [11] all cells infected at time tdie exactly at time t, thus in the one-step experiment their model predicts a perfect step-like result, in disagreement with experimental data (Fig. 1). Thus we expect the model by Jones et al. [11] to yield faster speeds than our model for two reasons: (i) they neglect the di¤usive delay; and (ii) they neglect the fact that the death of some cells takes longer than after infection. Replacing Dby 6 Deff (as explained in Sec. III), the model by Jones et al. is (see Eqs. (2.2) in Ref. [11]) @[I](r; t) @t =k1[V](r; t)[B](r; t) k1[V](r; t )[B](r; t );(17) @[V](r; t) @t =Deff @2[V](r; t) @r2k1[V](r; t)[B](r; t) +Y k1[V](r; t )[B](r; t );(18) @[B](r; t) @t =k1[V](r; t)[B](r; t):(19) Note that this is the same as the model by Yin et al. [Eqs. (12)-(14)], with k2I(r; t)replaced by k1B(r; t )V(r; t ). By following again the same method as in Sec. II, we …nd that the characteristic equation for the model due to Jones et al. [11] is 2c +Yec 1= 0:(20) Note that, in fact, the equation for @[I](r;t) @t above is not necessary to compute this speed, since [I]does not appear in the other two equations of the model by Jones et al. [11]. In Fig. 2 (rhombuses) we have also included the predictions of the model by Jones et al., for the same parameter values used in our model. We see in Fig. 2 that their model [11] predicts faster speeds than our model, as expected. Moreover, they are faster than the experimental speeds. For the wild strain, our model is consistent with the experimental range. For the mutants p001 and p005, the mean speeds predicted by our model are also closer to the experimental means (although the error bars are larger for the model by Jones et al. [11], because the speed depends strongly on k1). There is one more time-delayed model of virus front spread, due to Gourley et al. [13]. It is very similar to that by Jones et al. [11], discussed above, but it assumes an additional, natural death process only for infected cells (with rate Iand unrelated to virus infection) that decreases the number density of infected cells after time by a factor eI[13]. Although no biological reason was given in Ref. [13] why an additional death process might a¤ect only the infected cells (and not the uninfected ones), for completeness we next explore whether this model by Gourley et al. [13] changes the results of the model by Jones et al. [11] or not. Since this model by Gourley et al. [13] includes an additional death process for the infected cells, intuitively we expect that it could yield slower speeds than the model due to Jones et al. [11]. For the experimental conditions corresponding to the speeds that we analyze in the present paper (Fig. 2), the model proposed by Gourley et al. is (see Eqs. (1.1), (2.1) and (4.1) in Ref. [13]) @[I](r; t) @t =k1[V](r; t)[B](r; t) eIk1[V](r; t )[B](r; t );(21) @[V](r; t) @t =Deff @2[V](r; t) @r2k1[V](r; t)[B](r; t) +YeIk1[V](r; t )[B](r; t );(22) @[B](r; t) @t =k1[V](r; t)[B](r; t);(23) where we have neglected cell reproduction because in the experiments we want to explain, the cells were in the stationary growth phase before the arrival of viruses (as explained in Sec. II). We have also neglected virus death because it is negligible [16]. We do not include di¤usion of uninfected or infected cells because bacteria are immobilized in agar in these experiments (as mentioned in Sec. II). By following again the same method, the characteristic equation in the model due to Gourley et al. is 2c +Yec eI1= 0:(24) Again, in fact the equation for @[I](r;t) @t above is not necessary to compute this speed, since [I]does not appear in the other two equations of the set. In Fig. 2 (plotted as white squares) we have also included the predictions of this model by Gourley et al. [13] using the experimental value I= 0:4h1(from Fig. 7 in Ref. [17]). It is seen that its predictions are slower (as expected) but almost the same as those of the model by Jones et al. [11]. The speeds from both models are faster than the experimental ones. Finally, it is worth to note that, in situations where infected cells exit that class due to some other form of interaction, it would be necessary to modify our model. For example, for an additional, natural death process with exponential dynamics for infected cells, the righthand side in Eq. (4) would include an additional term I[I] (r; t)and Eq. (5) should be modi…ed accordingly. VII. CONCLUSIONS We have proposed a new reaction-di¤usion model with an eclipse time, that satisfactorily explains the experimental results of T7 virus plaques on E. coli. This improvement over previous models has been attained by means of the careful modi…cation of one of the evolution equations, which lacked biological signi…cance. Indeed, some previous models [7–9] assumed that the death rate of infected cells is proportional not only to their density, but also to the free space [Eq. (2)], which is not biologically reasonable. In contrast, the new model assumes that the death rate is proportional only to the density of infected cells (which begin to die after a time lag , corresponding to the eclipse phase of Fig. 1) [Eq. (3)]. Thus our new model is more reasonable biologically. Moreover, our new model agrees reasonably well with experimental data, in contrast to the classical model without delay or eclipse time due to Yin and et al. [4, 5]. 7 It is important to stress that Yin and co-workers already noted that their model was too fast for realistic parameter values, and only by …tting three parameters could it yield su¢ ciently slow speeds to agree with the experimental ones (Fig. 3 in Ref. [4]). In contrast, here we have not …tted any parameter but used realistic values, i.e. all parameter values we have applied have been obtained from independent experiments. Other authors took into account the role of the eclipse or delay time , but only in the reactive and not in the di¤usive process [11, 13], and they assumed the same eclipse time for all viruses. Those models yield faster speeds than the experimental ones. Also, we stress the importance of using realistic terms to modelize the interactions, e.g. the death process of infected cells (i.e. the release of viral progeny). Since the propagation of viruses is an active …eld of study in biophysics and medicine, having an underlying theory that is both mathematically and biologically sound is of special relevance. Furthermore, we have found that the results agree with experiments. By means of the detailed analysis of a simple mathematical model, we have aimed to demonstrate that such physical models are able to explain the spatial dynamics of virus infections. Certainly, in order to have a more comprehensive understanding of the problem, extensive data gathering for several viruses and environments should be undertaken. Appendix A: Time-delayed di¤usion In order to make this paper as self-contained as possible, here we include a brief derivation of Eq. (5). The derivation below (see Refs. [12, 18] for details) was originally proposed for human populations [12] and later applied to viruses [7–9]. During a time interval equal to the eclipse time (estimated from Fig. 1 in our case), the virus concentration changes both due to the reactive processes (1) and to dispersal. We …rst calculate the former change by using a Taylor series, [[V] (x; y; t +)[V] (x; y; t)]r=@[V] @t r +2 2 @2[V] @t2r + =F +2 2 @[V] @t r + (25) where the subindex rdenotes reactive processes, and F([V]) @[V] @t ris given by Eq. (7) according to the corresponding experiments (see Sec. II and Ref. [7]). Secondly, the change due to dispersal can be calculated by de…ning the dispersal kernel (x;y)as the probability per unit area that a virus initially placed at (x+ x; y + y)has moved to (x; y)after a time interval . Thus, [[V] (x; y; t +)[V] (x; y; t)]d =ZZ [V] (x+ x; y + y; t)(x;y) dxdy [V] (x; y; t):(26) In a system involving both reactive and dispersal processes, we add up their contributions [V] (x; y; t +)[V] (x; y; t) = ZZ [V] (x+ x; y + y; t)(x;y) dxdy [V] (x; y; t) + F +2 2 @F @t r + (27) Assuming that the kernel is isotropic, i.e., (x;y) = (), with  = q2 x+ 2 y;and Taylor-expanding Eq. (27) up to second order in time and space, @[V] @t + 2 @2[V] @t2=D@2[V] @x2+@2[V] @y2+F+ 2 @F @t r ; (28) where D=h2i 4=h2 xi 2=h2 yi 2is the di¤usion coe¢ - cient. For large distances r=px2+y2from the inoculation point of viruses ((x; y) = (0;0)), @2[V] @x2+@2[V] @y2'@2[V] @r2 and Eq. (28) is the same as Eq. (5), with Fgiven by Eq. (7) and Dreplaced by Deff (the reason for the latter change is explained in Sec. III). Thus the terms proportional to in Eq. (5) arise simply from a secondorder Taylor expansion. If the role of the eclipse time is neglected ('0), Eq. (28) reduces to the non-delayed or classical model used by Yin and co-workers [4, 5], namely (see Eq. 13) @[V] @t =D@2[V] @x2+@2[V] @y2+F: (29) In general, adding up the reactive and di¤usive contributions (as done in Eq. (27)) may not be exact [19–24] and this point is taken into account by the so-called sequential or cohabitation models (see especially Ref. [19], Fig. 1 of Ref. [20] and Fig. 17 of Ref. [21]). However, for virus infections cohabitation models yield almost the same results as non-cohabitation (or additive) models [23]. Thus in the present paper, we do not take the cohabitation e¤ect into account for mathematical simplicity (the predicted speeds in Fig. 2 would be the same, so there is no need to use more complicated equations). Let us mention that, in contrast to virus infections, for human waves of advance the cohabitation e¤ect is not negligible (and a more important e¤ect still is due to the shape of dispersal kernels) [20, 25]. Such more precise models lead to the ballistic speed for fast reproduction [22, 26], as they should [22, 24, 26]. However, for virus 8 infections those corrections are not necessary. In conclusion, the reaction-di¤usion Eq. (5) has the microscopic derivation above and recent criticisms [24] are irrelevant. For '0and F= 0;this also provides a valid derivation of Fickian di¤usion (Eq. (29) for F= 0). It is very important to stress that mathematical arguments [24] are not enough to establish whether a given equation is valid or not, because this depends on the system considered, and must thus be checked by using reactive functions, parameter values and initial conditions appropriate to the experimental setup (for example, to describe the growth of virus plaques, a model with only a pure death process is not realistic, and therefore irrelevant). As another example of this, reaction-di¤usion with Fickian di¤usion [Eq. (29)] can be applied if the delay time is negligible, which may be justi…ed for some biological species but not for viruses. This is clearly seen in Fig. 2, by comparing our model to the classical or non-delayed one (12)-(14) used by Yin et al. [4, 5], which is based on Eq. (29). At the other extreme, the second-order approximation would obviously fail for large [7, 24, 27] and, if this happened, additional terms in the Taylor expansions above would be necessary. However, this is not our case (Appendix B). Finally, Fickian di¤usion [Eq. (29) with F= 0] can be applied if the di¤usive delay time is su¢ ciently small, and is useful in many situations (not in our case). Thus parameter values must be examined to choose the appropriate equation for each experiment. Mathematical arguments are not enough, because an equation may be useful to describe some experiments but not others. For completeness, in Appendix B we extend the derivation above to in…nite order and …nd that the results are similar to those above and in the main paper (second order). Appendix B: Full time-delayed equation As shown in Appendix A, Eq. (5) is in fact an approximation, because it includes only terms up to second order from the Taylor expansions. The virus density [V]rapidly changes on a scale of time smaller than 15 min (because the increases in Fig. 1 take 6 minutes or less). This could therefore lead to errors in the front speeds obtained in Secs. IV-V. In this appendix we prove that this is not a problem by considering the full time-delayed equation (see Ref. [27], Eqs. (16) and (21)), 1 X n=1 n n! @n[V](r; t) @tn= 1 X n=1 (2Deff )n (2n)! @2n[V](r; t) @r2n + 1 X n=1 n n! @n1F(r; t) @tn1r ;(30) instead of its approximation, Eq. (5), together with Eq. (6) and our new Eq. (4). Then, repeating the same steps as in Sec. II we get the following characteristic equation (which replaces Eq. (8)), ec cosh p2e + 1c + ec  =Y c +1ec :(31) Repeating the calculations leading to Fig. 2, but using Eq. (31), we obtain that the di¤erences are very small. Indeed, the error between the second-order approximation and full time-delay equation is lower than 3% for the three strains of the T7 virus. Thus, the use of the second-order approximation in Secs. II-V is valid. Acknowledgments This work was partially funded by ICREA (Academia award) and the MINECO (projects SimulPast-CSD201000034, FIS-2009-13050 and FIS-2012-31307). [1] Weinbauer M. G., FEMS Microbiol. Rev.,28 (2004) 127. [2] Wodarz D., Hofacre A., Lau J. W., Sun Z., Fan H., et al., PLoS Comput. Biol.,8(2012) 1002547. [3] Koch A. L., J. Theor. Biol.,6(1964) 413. [4] Yin J., and McCaskill J. S., Biophys. J., 06 (1992) 1540. [5] You L., and Yin J., J. Theor. Biol., 200 (1999) 365. [6] Fort J., J. Theor. Biol.,214 (2002) 515. [7] Fort J., and Méndez V., Phys. Rev. Lett., 89 (2002) 178101. [8] Ortega-Cejas V. , Fort J., Méndez V., and Campos D., Phys. Rev. E,69 (2004) 031909. [9] Amor D.R., and Fort J., Phys. Rev. E, 82 (2010) 061905. [10] Yin J., J. Bacteriol., 175 (1993) 1272. [11] Jones D. A., Smith H. L., Thieme H. R., and Röst G., SIAM J. Appl. Math.,72 (2012) 670. [12] Fort J., and Méndez V., Phys. Rev. Lett.,82 (1999) 867. [13] Gourley S. A., and Kuang Y., SIAM J. Appl. Math.,65 (2005) 550. [14] Ebert U., and van Saarloos W.,Physica D,146 (2000) 1. [15] Yin J., Amn. N. Y. Acad. Sci.,745 (1994) 399. [16] De Paepe M., Taddei F., PLoS Biology,4(2006) e193. [17] Zobell C. E., and Cobet A. B., J. Bacteriol.,84 (1962) 1228. [18] Isern N., and Fort J., Phys. Rev. E,80 (2009) 057103. [19] Fort J., Pérez-Losada, J., and Isern, N. Phys. Rev. E,76 (2007) 031913. [20] Isern N., Fort, J., Pérez-Losada, J. J. Stat. Mech.: Theor. & Exp., 2008 (2008) P10012. [21] Fort J., Pujol, T., Rep. Progr. Phys.71 (2008) 086001. [22] Fort J., PNAS, 109 (2012) 18669. [23] Amor D. R., and Fort J., Physica A, 416 (2014) 611. [24] Méndez, V., Campos, D., and Horsthemke, W., Phys. 9 Rev. E, 90 (2014) 042114. [25] Fort J., J. R. Soc. Interface, 12 (2015) 20150166. [26] Fort J., Pérez-Losada J., Suñol J. J., Escoda, L., and Massaneda, J. New J. Phys., 10 (2010) 043045. [27] Fort J., and Méndez V., Phys. Rev. E,60 (1999) 5894.