Fronts with Continuous Waiting-time Distributions: Theory and Application to Virus Infections

We generalize to arbitrary waiting-time distributions some results which were previously derived for discrete distributions. We show that for any two waiting-time distributions with the same mean delay time, that with higher dispersion will lead to a faster front. Experimental data on the speed of virus infections in a plaque are correctly explained by the theoretical predictions using a Gaussian delay-time distribution, which is more realistic for this system than the Dirac delta distribution considered previously ͓J.


I. INTRODUCTION
The speed of reaction-diffusion fronts is a very important problem in many areas of physics, such as flame propagation, solidification, superconductors, and biophysics ͑for some recent reviews, see ͓1,2͔͒.On the other hand, timedelayed transport is observed in many different situations ͓3,4͔.For instance, time-delayed diffusion-reaction fronts are observed in many systems ͑specially, biophysical ones͒ because the individuals ͑or particles͒ which diffuse and reproduce ͑or react͒ have a non-negligible waiting time between successive jumps ͓1,2,5-11͔.In the first applications of this topic, the approximation was made that all individuals have exactly the same waiting time, i.e., that the waiting-time distribution is a Dirac delta ͓5,6,8͔.But of course, a Dirac delta is a very specific, highly idealized description, and should not be expected to be valid in general.For example, for the Neolithic transition population fronts, the available data were recorded as a discrete distribution of several possible values of the waiting time ͑and not a single one͒.For this reason, in Ref. ͓7͔ we extended the time-delayed theory of the Neolithic transition ͓5͔ to the case of several possible, discrete waiting times.In Ref. ͓7͔, we showed that there is an effective waiting time T ˜, and that this is a very important parameter to compute and understand front speeds.However, in other systems a discrete distribution of delays ͑assumed in Ref. ͓7͔͒ will not be realistic.As we shall show in the present paper, this happens for virus infection fronts: The waiting time distribution is not a sum of Dirac deltas ͑discrete distri-bution͒, but a Gaussian one ͑i.e., a continuous distribution͒.Therefore in order to deal correctly with virus infection fronts, it is necessary to extend the approach in Ref. ͓7͔ to continuous waiting-time distributions.Indeed, this is the main experimental motivation for the present work.But as explained above, from a theoretical perspective it is also important to build a framework which makes it possible to compute the effective waiting time T ˜and front speed c, and is applicable to all waiting time distributions ͑and not only to discrete ones, which was the case considered in Ref. ͓7͔͒.This problem is tackled in Sec.II and is applied to virus infections in Sec.III.Finally, Sec.IV is devoted to concluding remarks.

II. THEORY
For the sake of clarity, we will refer to bacteriophages in this section.This will make it easier to apply our results and compare to experimental data in the next section.Bacteriophages are viruses that infect bacteria.Some time after infection, the bacteria burst and release new virus progeny.However, the framework presented in this section is applicable to any kind of individuals or particles performing a random walk with a distribution of waiting times and undergoing biological reproduction/replication or chemical reactions.For example, in addition to bacteriophages, the same model can be applied to vesicular stomatitis viruses ͑VSV͒, which replicate on mammalian or insect cells ͑not on bacte-ria͒.In VSV, the cells do not burst, but there is again a waiting time distribution due to the difference in release times of viruses from the infected cells.
In the bacteriophage plaque experiments, viruses diffuse and reproduce in a two-dimensional medium.Virus diffusion takes place in "jumps" from cell to cell.Once a virus V reaches a cell ͑i.e., a bacterium B͒, it is adsorbed and infects it ͑this is represented by the reaction V + B → I͒.The virus then reproduces inside the infected bacterium I. Some "waiting" time after the arrival of the virus, the infected bacterium I bursts and a number Y of new viruses ͑progeny͒ leave it ͑reaction I → YV, where Y is called the yield͒.Bacteria are immobilized and, therefore, do not diffuse ͑see ͓8͔ and references therein͒.Let V͑x , y , t͒ stand for the virus number density per unit area in position ͑x , y͒ at time t.Assume that initially ͑t =0͒ we have V͑x , y , t͒ = V 0 ␦͑x͒␦͑y͒, i.e., that initially a concentration V 0 of viruses is inoculated at a point which is taken as the origin ͑this is indeed how the experiments considered in Sec.III are performed in the laboratory͒.Let ͑⌬x , ⌬y͒ stand for the probability that a virus moves with coordinate lengths −⌬x, −⌬y when performing a jump for a cell to another one.Let ͑T͒ stand for the probability that from the time a virus reaches a cell, a waiting time T elapses until its progeny leave the cell.Let F͑x , y , t͒ stand for the number of new viruses that appear per unit time and area, due to the adsorption and replication reactions above ͑i.e., V + B → I → YV͒ ͓8͔.If ds P V ͑x , y , t͒ stands for the number of viruses per unit area that reach an area ds centered at ͑x , y͒ at time t, we have ͓7͔ The density V͑x , y , t͒ of viruses per unit area centered at x͑x , y͒ at time t is clearly given by the viruses that have arrived at ͑x , y͒ at some earlier time and still not left, namely ͓7͔ ͑x,y,t͒ = ͵ 0 t dtЈP͑x,y,tЈ͒⌿͑t − tЈ͒, ͑2͒ where ⌿͑t − tЈ͒ is the probability that any particle rests for at least a time interval t − tЈ before performing the next jump, obviously ͓7͔ It is straightforward to repeat the steps leading to Eq. ͑11͒ in Ref. ͓7͔ and see that here the same general equation holds, i.e., ͓12͔

͑4͒
where V ˆ͑k x , k y , s͒ and F ˆ͑k x , k y , s͒ are the Fourier-Laplace transforms of V͑x , y , t͒ and F͑x , y , t͒, ˆ͑s͒ is the Laplace transform of ͑T͒, and ˆ͑k x , k y ͒ is the Fourier transform of ͑⌬x , ⌬y͒.We now proceed as follows.
͑i͒ As in Ref. ͓7͔, we assume that the space kernel is isotropic, i.e., ͑−⌬x , ⌬y͒ = ͑⌬x , ⌬y͒ = ͑⌬x ,−⌬y͒ = ͑⌬y , ⌬x͒, which using the normalization of probability ͓͐ −ϱ ϱ d⌬x͐ −ϱ ϱ d⌬y ͑⌬x , ⌬y͒ =1͔ leads to where O͑⌬ ជ 3 ͒ stands for terms of third and higher powers of ⌬x and ⌬y.This approximation will be valid assuming that the dispersal kernel ͑⌬x , ⌬y͒ is appreciably different from zero only for sufficiently small jumps ͑⌬x Ӎ 0, ⌬y Ӎ 0͒.Otherwise, the second-order or "diffusion" approximation above would break down ͑leading to what is called long-range dispersal in ecology ͓13-15͔͒.͑ii͒ In Ref. ͓7͔ we assumed that the waiting time probability distribution was discrete, i.e., a sum of Dirac deltas: That assumption was motivated by the available experimental data of ͑T͒ in the application considered there ͑namely, the Neolithic transition͒.Here, we could be tempted to use a Gaussian distribution ͑which is motivated by the experiments that we will consider in Sec.III͒.However, we will now show that we can derive very interesting results without assuming any specific distribution.Instead, and analogously to Eq. ͑5͒, we simply assume that the waitingtime probability distribution ͑T͒ is appreciably different from zero only for sufficiently small values of the waiting time T, so that we can use again a second-order Taylor expansion as follows: Combining the three previous equations up to second order, we obtain that where we have defined T ˜can be called the effective delay time.It was first introduced in Ref. ͓7͔, but for a very special and simple distribution function ͑T͒ ͑namely, a sum of Dirac deltas͒.We now see that T ˜is a very important parameter in general.The meaning of the effective delay T ˜can be understood by rewriting Eq. ͑8͒ as where is the dispersion of the waiting-time distribution.Thus for any two waiting-time distributions with the same mean delay time ͗T͘, that with higher dispersion will have a lower effective delay time T ˜and, therefore, a faster front.In other words, some viruses jump sooner ͑due to the higher disper-sion͒ and make the infection front move faster.Conceptually, this effect is somehow similar to long-range dispersal in ecology ͓13-15͔.There, a few seeds dispersing large distances can lead to a much faster front.Here, a few viruses dispersing sooner can also lead to a faster front.

͑11͒
Equations ͑7͒ and ͑11͒ are analogous to Eqs. ͑19͒ and ͑20͒ in Ref. ͓7͔, but have been here derived for an arbitrary distribution function ͑T͒.It is very interesting that reactiondiffusion systems follow an HRD equation, not only for the particular distribution of waiting times considered in Ref. ͓7͔ ͑a sum of Dirac deltas͒, but for any general waiting-time distribution ͑T͒, as shown above.
Just to summarize: originally ͓5͔, the validity of an HRD equation up to second order in waiting-time random walks was shown under the assumption of a single possible value T 1 for the waiting time, i.e., for ͑T͒ = ␦͑T − T 1 ͒ ͑with the result that T ˜= T 1 , using the notation in the present paper͒.Later on ͓7͔, an HRD equation was derived for a discrete set of waiting times, i.e., for ͑T͒ = ͚ i=1 N p i ␦͑T − T i ͒.Now, an HRD equation has been derived for any possible waiting time distribution ͑T͒ ͓with the effective delay T ˜appearing in the HRD Eq. ͑11͒ given by Eq. ͑8͔͒.

A. General model
As discussed in Ref. ͓8͔, for virus infections the diffusion coefficient D must be replaced by an effective one, D ef f , to take into account the presence of bacteria ͑which hinder virus diffusion͒ as where f = B 0 / B max is the concentration of bacteria relative to its maximum possible value ͑B 0 is the initial bacteria concentration far from the inoculation origin, and it depends on the initial nutrient concentration͒ ͓8͔.The parameter x takes proper care of the bacterial shape ͓8͔.Also, the evolution equations for the uninfected bacteria number density B͑x , y , t͒ and infected bacteria density I͑x , y , t͒ are ͓8͔ where k 1 is the rate constant of the virus adsorption reaction ͓V + B → I͔, k 2 the rate constant of the infected bacteria lysis reaction ͓I → YV͔, and I max the saturation density of infected cells.Finally, the reactive term for viruses is ͓8͔ Therefore using the HRD Eq. ͑11͒ we have a system of three simultaneous partial differential equations, which generalize those considered in Ref. ͓8͔.In Eq. ͑16͒, instead of the mean waiting time ͑used in Ref. ͓8͔͒, now the effective delay T ˜appears, and it must be computed using Eq.͑8͒.This takes care of the whole waiting-time distribution function ͑which was not done in Ref. ͓8͔͒.The solution obtained by linearization in the front frame by requiring the determinant of the matrix corresponding to the linearized form of Eqs.͑16͒ to vanish.This solution will be the same as that in Ref. ͓8͔ with the mean waiting time replaced by the effective one T ˜, namely ͓see Eq. ͑8͒ in ͓8͔͔ where 1 ϵ k 1 B 0 / k 2 and T ˜ϵ T ˜k2 are dimensionless parameters.This equation can be solved numerically in order to find out the dimensionless front speed c ¯ϵ c / ͱ D ef f k 2 such that c ¯= min Ͼ0 ͓c ¯͔͑͒, where c ¯͑͒ is given by characteristic equation ͑17͒.But before applying this equation we must estimate T ˜from Eq. ͑8͒, derived in the previous section.In order to do so, we first have to determine the waiting-time distribution function.This was not done in Ref. ͓8͔ because in that very simple model, all viruses were assumed to have exactly the same waiting time, i.e., we approximated the waiting-time distribution to a Dirac delta centered at the mean waiting time of the viruses.

B. Bacteriophage T7 infecting E. Coli
In order to analyze the waiting-time distribution for a real system, in Fig. 1 we reproduce the so-called one-step growth of the virus T7 infecting E. Coli bacteria ͑see also Fig. 1 in Ref. ͓8͔ and its caption͒.This experiment refers to an homogeneous medium of cells infected at t = 0.If all viruses took exactly the same time to kill a cell and reproduce, Fig. 1 would be a step function, and the waiting time distribution would be a Dirac delta.Instead, the gradual rise in the virus concentration in Fig. 1 indicates that it takes a different time for each virus to kill the cell it has infected and reproduce.This interpretation of one-step curves is well-known in virology ͓16͔.For the case in Fig. 1, we see that the range of waiting times ͑i.e., the rise in the curve͒ is between 14 and 23 min, approximately.Instead, Ref. ͓8͔ assumed that the waiting time of all viruses is exactly equal to the mean value in Fig. 1 ͑i.e., 18.4 min͒.Intuitively, it is clear that there is no reason to expect that this approximation will lead to realistic predictions ͑the width of the rise in Fig. 1 2 are always higher than 13% ͑and even higher than 65% for sufficiently small values of f͒.Therefore viruses that kill an infected cell sooner than average ͑T = 14 min͒ tend to make the infection front move substantially faster than would be expected from the single, mean waitingtime model in Ref. ͓8͔ ͑T = 18.4 min͒.But on the other hand, computing the time derivative of the main Fig. 1 ͑inset in Fig. 1͒ we note that very few viruses kill the cell in only 14 min.Therefore in order to make trustable, quantitative predictions we need to take into account the detailed shape of the waiting time distribution.We will then be able to predict the virus front speed using the general model presented in Sec.II.
As mentioned above, the experimental data in Fig. 1 were obtained for a homogeneous medium of cells infected at t = 0.Then, since each virus disappears and gives rise to a progeny of Y viruses after a time T with probability ͑T͒, obviously the concentration of viruses will in that experiment evolve according to ͓17͔ so that the waiting-time probability distribution can be obtained from the curve in Fig. 1 as The inset in Fig. 1 shows the time derivative of the main Fig. 1 ͑full curve͒, and a Gaussian fitted by least-squares ͑dotted curve͒.It is seen that a Gaussian is a very good description of the waiting-time distribution of these viruses ͑in contrast, the approximation in Ref. ͓8͔ considered a vanishing distribution for all values of T except 18.4 min, i.e., a Dirac delta distribution͒.Therefore here we will use a Gaussian waitingtime distribution, so that the normalization constant ͓i.e., the value of A such that ͐ 0 ϱ dT͑T͒ =1͔ and the mean squared waiting time ͗T 2 ͘ are, respectively, where Erf͓z͔ϵ͑2/ ͱ ͒͐ 0 ϱ exp͓−t 2 ͔dt is the error function.The former results become much simpler if we consider the special case that all viruses have a waiting time substantially different from zero.In other words, if we consider the case in which the time between the arrival of a virus and the departure of its progeny is not negligible for any of the viruses ͑below we shall see that this is indeed realistic͒.Intu-FIG.1. Virus concentration ͓in Formation-of-Plaque Units ͑FPU͒ per ml͔ vs time in a homogeneous medium of cells infected at t = 0.The fit to the main plot is a logistic ͑see main Fig. 1 in Ref. ͓8͔͒.Here we use it again because its time derivative ͑inset, full curve͒ makes it possible to note that a Gaussian ͑dotted curve͒ is a good description to the waiting-time distribution of the T7 bacteriophage.
itively, we may express this condition by means of the mathematical inequality

͑24͒
In this special case, the Gaussian function in Eq. ͑20͒ is approximately zero for T Ͻ 0, and we may approximate the normalization condition as follows: and Eq.͑22͒ becomes the very simple expression which below we shall see that is realistic and very useful.The Gaussian curve fitted to Eq. ͑20͒ is shown as a dotted curve in the inset in Fig. 1.It has the parameter values ͗T͘ = 18.38 min and B = 1.634 min.Using these values into Eq.͑22͒ yields ͗T 2 ͘ = 339.1 min 2 .The same result can be found from the approximation ͑27͒ ͓because for these values exp ϫ͓−͑͗T͘ / B͒ 2 ͔ϳ10 −55 , so that the condition ͑24͒ holds͔.Then, using Eq.͑8͒ we can estimate the effective waiting time, T ˜= 18.31 min, which is very similar to the value of the mean waiting time ͗T͘ = 18.38 min obtained above ͑from Fig. 1, inset͒.Using the rest of the parameter values for our sys-tem from Ref. ͓8͔, we obtain, solving Eq. ͑17͒ numerically, the virus front speed predictions shown in Fig. 3 ͑curves͒, which agree well with the experimental data ͑symbols͒.For this system, the predictions are very similar to those from the single-delayed model ͓Fig.3͑a͒ in Ref. ͓8͔, lower curves͔, and we have checked that the same happens for B max =10 8 ml −1 ͓Fig.3͑b͒ in Ref. ͓8͔͔.This happens because we have obtained that T ˜Ӎ͗T͘ for this system.However, the arguments above Eq.͑18͒ and Fig. 2 show clearly that it was necessary to perform this analysis.We can now conclude that physical models can explain the virus front experiments, contrary to the widespread misbelief that they are driven by unknown biological factors ͑see ͓8͔ and references therein͒.
The simplest model is that in Ref. ͓8͔, where a single waiting time was assumed.Here, we have seen that taking into account the shape of the waiting time distribution yields similar results.Before leaving this example, let us stress again that the former section contains a general framework that can be used for any diffusive system with a distribution of waiting times ͓18͔.

C. Influence of the waiting-time distribution on the front speed
Up to now we have applied our new theoretical results to the experimental data in Ref. ͓8͔.This analysis was necessary in order to have trustable predictions because a Gaussian distribution is realistic in this case ͑Fig. 1, inset͒.We have found that T ˜Ӎ͗T͘ for the experimental data in Ref.Both curves are based on that model, but the dotted one uses the lower value of the range of waiting times implied by Fig. 1, i.e., T = 14 min, whereas the full curve uses the mean value T = 18.4 min.Differences between both curves are substantial, showing the need to build models that take care of the shape of the waiting time distribution.This is done in the present paper.measurements for such a case, it is interesting to analyze how much would the Gaussian distribution ͑20͒ have to change from a Dirac delta so that the front speed would be altered substantially.This is the question we tackle in this section.
Provided that the condition ͑24͒ holds, Eq. ͑27͒ shows that the effective and mean waiting times ͑T ˜and ͗T͘, respec-tively͒ are related to the width B of the Gaussian distribution ͑20͒ as so that the dispersion of the waiting-time distribution ͑10͒ can be written as From this and Eq.͑9͒, we can estimate the relative difference between the effective delay T ˜and the mean delay ͗T͘ as As we shall now see, this equation is very useful to predict whether the results of the Gaussian and the Dirac-delta distributions will be substantially different or not.In passing, let us just recall that, as it is clear from the previous equation, the effective delay T ˜is always smaller than the mean delay ͗T͘, and the physical meaning of this mathematical inequality is that the viruses ͑or particles͒ with low values of waiting time T ͑the left part of the Gaussian distribution͒ tend to drive the front faster than if only viruses with the mean delay ͗T͘ were considered ͑Dirac distribution͒.
From Eq. ͑30͒ we can predict the difference between the effective delay T ˜and the mean delay ͗T͘.Let us consider three important examples.
͑i͒ The case in Ref. ͓8͔ ͑a Dirac delta waiting-time distri-bution͒ is recovered here as a special case, namely that of a Gaussian distribution with vanishing width ͑B =0͒, which from Eq. ͑30͒ yields T ˜= ͗T͘, so that only in this limiting, very special case does the HRD equation ͑11͒ derived in this paper reduce to that in Ref. ͓8͔.͑ii͒ For the case in the previous section ͑͗T͘ = 18.38 min and B = 1.634 min͒, Eq. ͑30͒ yields a difference between T ãnd ͗T͘ of only 0.4%.This is the reason for the similar predictions ͑see the former section͒ from the Dirac-delta model in Ref. ͓8͔ ͑for which ͗T͘ appeared in the HRD equation͒ and the Gaussian model here ͓for which T ˜appears in the HRD Eq. ͑11͔͒.
͑iii͒ Finally we can answer the question asked at the beginning of this section: How much does a Gaussian distribution have to differ from a Dirac delta to yield substantially different front speeds?The answer is given by the new Eq.͑30͒.For instance, consider a hypothetical Gaussian waitingtime distribution with B = ͗T͘ / ͱ 2. First of all, we easily check that the condition ͑24͒ holds approximately ͓20͔.Thus we can use Eq.͑30͒.In this way we obtain that the difference between T ˜and ͗T͘ is of 25% ͑thus B = 13 min and T = 13.79 min, if we use the same typical value ͗T͘ = 18.38 min as above͒.A difference of 25% between the relevant delay times in both models ͑namely, T ˜and ͗T͒͘ can be expected to yield substantially different front speeds.Indeed, we shall now give the results from both models for the typical values ͓8͔ f = 0.2, k 1 = 1.88ϫ 10 −9 ml/ min, B max =10 8 ml −1 , and the same other values of the parameters as above.First, the Dirac-delta model in Ref. Ref. ͓8͔ can be an oversimplified description that can lead to substantial errors.Note, once again, that the Gaussian model leads to a faster front because some viruses leave the infected cell sooner ͓note also that in example ͑ii͒ above, the Gaussian distribution had a much lower width B, leading to a substantially weaker effect͔.For some other typical values of f, B max , etc., the difference between both models increases but, as mentioned, there are no experimental data of front speeds available, so there is no point to present additional hypothetical results at this stage.Indeed, the previous example answers very clearly the relevant physical question, i.e., how the distribution has to depart from a Dirac delta to yield substantially different predictions.The important point is that Eq. ͑30͒, first derived in the present paper, provides a simple way to predict to what extent the predictions of the Diracdelta and the Gaussian models will differ.This equation shows that the predictions from both models will differ more the higher the Gaussian width B and/or the lower the mean waiting time ͗T͘, as was to be expected intuitively, because in both cases more viruses ͑or particles͒ will obviously jump sooner.tion to the time rate of change of the virus concentration in the case of nonhomogeneous media.One way to see this is that in such an instance, there are additional viruses that reach the point considered after t =0 ͑because of diffusion͒ and they will also replicate, but their progeny are not included in this equation.A second way to see the inapplicability of these equations in the plaque growth experiments ͑corresponding to Fig. 3 in the present paper͒ is that then V t=0 = 0 outside the small region initially occupied by viruses.This is why we use a logistic source term for the lysis reaction I → YV, namely the last term in Eq. ͑15͒ ͑note that it corresponds to the full line in the main Fig. 1, and it describes well the corresponding process, as explained in detail in Ref. ͓8͔͒.
is about 7 min, i.e., almost 40% of the mean value of 18.4 min, and is therefore not negligible͒.It is also interesting to use the simple, single waiting-time framework in Ref. ͓8͔ to justify that the distribution of waiting times should not be neglected a priori.Figure 2 presents an example of the predictions of the single waiting-time model in Ref. ͓8͔, but using the lower value implied by Fig. 1, i.e., T = 14 min ͑dotted curve͒ in addition to the mean value T = 18.4 min ͓full curve in Fig. 2, already plotted in Ref. ͓8͔, Fig. 3͑b͔͒.Differences in the front speed c in Fig.
FIG. 3. Predictions of the model in the present paper ͑curves͒ vs experimental data ͑symbols͒ for the front speeds of T7 viruses infecting E. Coli bacteria.
͓8͔and, therefore, that the predictions of the new model here are similar to those of the simpler model in Ref. ͓8͔.That simpler model assumed the same waiting time for all viruses ͑i.e., a Dirac delta distribution͒, whereas here we have considered a Gaussian distribution ͑Fig. 1, inset͒.As mentioned, for the experimental data in Ref. ͓8͔ the predictions of both models ͑the Dirac-delta model and the Gaussian one͒ are similar.But for other experiments, they may be substantially different.Although we are not aware of virus front speed FIG. 2. This figure presents the limitations of the single waiting-time model in Ref. ͓8͔.