Traveling fronts in self-replicating persistent random walks with multiple internal states

Self-activation coupled to a transport mechanism results in traveling waves that describe polymerization reactions, forest fires, tumor growth, and even the spread of epidemics. Diffusion is a simple and commonly used model of particle transport. Many physical and biological systems are, however, better described by persistent random walks that switch between multiple states of ballistic motion. So far, traveling fronts in persistent random walk models have only been analyzed in special, simplified cases. Here, we formulate the general model of reaction-transport processes in such systems and show how to compute the expansion velocity for arbitrary number of states. For the two-state model, we obtain a closed-form expression for the velocity and report how it is affected by different transport and replication parameters. We also show that nonzero death rates result in a discontinuous transition from quiescence to propagation. We compare our results to a recent observation of a discontinuous onset of propagation in microtubule asters and comment on the universal nature of the underlying mechanism.


I. INTRODUCTION
Transport and self-replication give rise to propagating fronts found across many systems in physics, chemistry, engineering, and ecology [1][2][3][4][5][6] . Examples range from opinion spreading 7 and forest fires 8 to tumor growth 9 and cell biological processes 10,11 . Diffusion is a convenient and widely-used model of the transport process 2,4,12 . This approach is often criticized because it neglects correlations or persistence between the motion at different times 13 . A less known, but perhaps more crucial, assumption behind the diffusion approximation is that the velocity of particles has no upper bound 14,15 . Indeed, the Gaussian distribution extends infinitely far from the starting position of a random walker. Such instantaneous transport could yield inaccurate or simply unphysical results in many situations of front propagation 1,15 . a) kishihara@pks.mpg.de b) korolev@bu.edu The simplest way to account for a finite transport velocity is to consider a particle that moves with a fixed velocity and reorients its direction at a constant rate. Such motion is known as persistent, correlated, or ballistic random walks and is described by the generalized telegraph equation in the continuum limit 13,[15][16][17][18] . Propagating fronts of self-replicating persistent random walks have been studied by a variety of methods including the Wentzel-Kramers-Brillouin approximation 19 , phase portrait analysis 4,15 , Hamilton-Jacobi approach 1 , and numerical simulations 20 . These studies showed that the diffusion approximation predicts front velocities that exceed the maximal velocity of the particle, while the telegraph equation correctly captures the upper bound on the front velocity 1,15 .
Most of the previous work focused on particles with just two states (left and right moving) and further assumed certain symmetries for reproduction, switching, and movement rates, which are not expected to hold in general 15,21,22 . Here, we develop a generalized framework for traveling fronts formed by persistent random walks with an arbitrary number of states. Switching between multiple states with different motility and replication is a realistic representation of the individual behavior of bacteria 23,24 , motile cells 25,26 , cancer cells 27 , and animals 15,28 . We show how to formulate an appropriate model, identify relevant parameter combinations, and compute the velocity of a traveling front. Similar to prior studies, our results are limited to the so-called "pulled waves" whose kinetics is controlled by the replication at the leading edge of the front 2, 29,30 .
For the two-state model, our method produces an explicit formula for the front velocity in terms of elementary functions. Unlike previous studies, we obtained the solution for the most general model, which contains ten distinct parameters. The exact solution greatly facilitates the exploration of the parameter space and clearly shows how different growth and transport states contribute to the overall motion of the front. We compare the exact solution to the diffusion approximation and pay special attention to important limiting cases. In particular, we derive conditions necessary for the onset of front propagation, determine the relative importance of different replication modes, and explain asymptotic behavior at small and large replication rates.
In many applications, the particles can not only replicate, but also die. Because this possibility has been rarely considered in the literature 31 , we thoroughly analyze the effect of death on front propagation. Our main finding is that the onset or cessation of propagation is accompanied by a discontinuous jump of the velocity. That is, as the replication rate increases, the velocity jumps from zero to a finite value when the population transitions from net death to net growth.
Velocity jumps have been recently observed in growing microtubule networks 32 . Similar to persistent random walks, microtubules and other biopolymers possess long-lived states in which the polymer either polymerizes or depolymerizes [33][34][35] . Many of such biopolymers also catalyze their own nucleation and thus self-assemble as traveling fronts 32,36,37 . We provide an intuitive explanation for velocity jumps in biopolymers and highlight the important differences between self-replicating particles and self-replicating polymers.

II. PHENOMENOLOGICAL DESCRIPTION: DIFFUSION EQUATION WITH ADVECTION AND GROWTH
To gain a conceptual understanding, we first describe a simple model of coupled growth and transport: a reaction-advection-diffusion equation, where n is the population density, D is the diffusion coefficient, U is the advection velocity, and g(n) is the per capita growth rate. The advection velocity is typically set to zero in the context of range expansions and biological invasions because there is no preferred direction for dispersal. However, in some situations, such as bacteria in the urinary tract or aquatic animals in rivers, advection velocity cannot be neglected because it accounts for the directed motion due to fluid flow in the environment 38,39 .
When g(0) > 0, any initial inoculation results in local growth until the density saturates atn such that g(n) = 0. In addition to local growth, the population spreads spatially producing two traveling fronts on each side of the population. These fronts translate in space with constant velocity V without changing their shape. In other words, n(t, x) depends only on x − V t rather than t and x separately. Here, we exclusively focus on pulled waves 2,40 , i.e expansions whose velocity is controlled only by the dynamics at the leading edge. Pulled expansions occurs when growth dynamics have negative or weakly positive feedback. For example, expansion is guaranteed to be pulled when g(n) is monotonically decreasing 2, 40 . The front velocities of pulled waves are obtained by linearizing equation (1) and are given by Each of the two fronts may move either leftward or rightward depending on the relative magnitude of the advection velocity U and the so-called Fisher velocity 2 g(0)D.
To fix ideas, let us consider logistic growth with death, which is commonly used in mathematical biology, Here, r is the growth rate at low population density, K is the carrying capacity, and d is the death rate. From equation (2), one can easily determine how front velocities depend on the model parameters and, for example, compute the critical growth rate r * = U 2 4D , which is necessary to propagate against the flow.
The onset of propagation can be either gradual when V changes continuously as model parameters are varied or sudden when V jumps from zero to a non-zero value. The gradual onset occurs only when d = 0, which is the case that has received most attention in the literature. A sudden onset is the generic case because it occurs whenever the death rate is positive. Indeed, imagine increasing the growth rate r from below d to above d. When r < d, the population size declines until the population becomes extinct. When r > d, the population grows, with the population translating with nonzero velocity U as soon as it is viable. In the present context, this observation is rather trivial, but, in more complex system such as biopolymers, it leads to important experimental signatures that can be used to elucidate biological function or mechanism 32 . This manuscript is concerned with developing the theory of the aforementioned phenomena in self-replicating persistent random walks (Sec. IV B) and polymers (Sec. V A).

III. A GENERAL FRAMEWORK FOR PERSISTENT RANDOM WALK OF REPLICATING PARTICLES WITH MULTIPLE STATES
As we discussed in the introduction, diffusion approximation permits arbitrary speeds of movement. Indeed, equation (2) predicts unbounded growth of V with g(0). This unphysical behavior can be eliminated by enforcing an upper bound on the particle velocity. Persistent random walk is the simplest and widely applicable models that constrains particle velocity. In this model, particles exist in multiple states each with a fixed velocity. The transitions between the states occur with constant, but potentially state-dependent rates.
Here, we formulate a general framework for the persistent random walk of replicating particles with an arbitrary number of states, and provide the exact solution for the front velocity. Throughout the paper, we use the term "states" to refer to particles that exhibit different velocity, switching, and replication rates. Depending on the specific application, "states" could be considered as species, phenotypes, particles, or behaviors.

A. Mechanistic formulation
Consider a one-dimensional system of replicating particles that can be in one of N distinct states. Under the assumptions stated above, the dynamics of this system is given by N coupled first-order differential equations: where n α is the density of particles in state α; v α is the particle velocity; d α is the death rate; f αβ ≥ 0 is the switching rate from state β = α into state α; and r αβ ≥ 0 is the replication rate at which particles is state β produce particles in state α. The switching rate out of state α is given by −f αα , which also equals γ =α f γα by the conservation of probability.
In general, all the rates can depend on all the population densities of the particles {n α }.
For pulled waves, the expansion velocities depend only on the values of these rate at low population densities (n α → 0), so we omit the density dependence in the following and implicitly assume that the values are taken in the limit of vanishing n α . In simulations, we need to impose density-dependence to prevent unbounded growth. We accomplished by setting r αβ (n) = r αβ (1 − n/K), where n = α n α is the total population density and K is the carrying capacity. All other rates had no density-dependence in our simulations.

B. Reduced formulation
The last three terms in equation (4) have the same functional form and can be combined for a more compact formulation where and δ αβ denotes Kronecker's delta. Note that the choice of Λ is somewhat constrained because the mechanistic rates: d α , r αβ , and f αβ need to be positive.
From the reduced formulation, it is clear that front velocities and other quantities will depend on the mechanistic rates only through Λ αβ . In many situations, it is however important to understand how changes in specific mechanistic rates affect population dynamics. Thus, we need a way to decompose Λ αβ into its constitutive parts. This decomposition is, however, not unique. When needed, we overcome this ambiguity by imposing additional requirements. Specifically, we require that β f αβ = 0, i.e. switching between states conserves the number of particles and set f αβ = Λ αβ for α = β. The diagonal terms g α then accounts for both death and replication, which is assumed to produce only particles of the same type as the parent. Thus, Λ αβ = g α δ αβ + f αβ , where g α can take take both positive and negative values.

C. Asymptotic solution for pulled fronts
Velocity of the traveling fronts is the key property of a population of self-replicating persistent random walks. Without loss of generality, we focus on the velocity of the right edge of the population. This velocity can be computed following the procedure described in Ref. 2 . The detailed derivation is provided in Appendix A 1, but the outline of the calculation is summarized below. As noted earlier, our solution applies to pulled fronts whose velocity can be computed by linearizing the dynamical equation and setting Λ αβ to Λ αβ (0).
We seek the long-time asymptotics for the traveling front, characterized by a constant velocity and shape. Following linearization, the dynamical equation is Fourier transformed in space and Laplace transformed in time. The Fourier variable is denoted as k and Laplace variable as s. The resulting system of linear algebraic equations is solved by the standard methods, and the inverse Fourier and Laplace transforms are performed. The resulting integrals are evaluated in the long time limit using saddlepoint approximation.
The key quantity that arises in the above-described calculation is the following matrix where parameter κ is real (the relevant Fourier mode has purely imaginary k, so we introduced κ = −ik > 0). The eigenvalues of this matrix, s (l) , describe available saddle points. The saddle point that predicts the largest velocity dominates the long-time dynamics.
For each eigenvalue, the corresponding front velocity is determined by solving with respect to κ. We denote the solution of equation (8) by κ (l) f . This quantity specifies the spatial decay rate of all particle densities in the comoving reference frame n α ∼ e −κ (l) f is known the front velocity is given by The actual velocity and spatial decay rate are given by l that predicts the largest V (l) : For a two-state system with velocities v 1 and v 2 , the front velocity of the right side of the population can be computed in closed form: where sgn is the sign function. The detailed derivation of this result is provided in Appendix A 2. We note that the above formula contains a lot of nontrivial information. Indeed, the model has six independent parameters and only one of them could be set to unity by rescaling time.
In the absence of growth and death, our result for V reduces to the average velocity (V =v), which is defined as followsv where f 1 denotes f 21 and f 2 denotes f 12 , i.e. f 1 and f 2 are the rates of switching out of states 1 and 2 respectively.

IV. PERSISTENT RANDOM WALK OF REPLICATING PARTICLES
To understand the implications of equation (11), it is convenient to recast the model into an intuitively interpretable form of persistent random walks; see figure 2a. We assume that random walks can move either rightward or leftward with speeds v + and v − . Switching between the two states occurs at rates f + ≥ 0 and f − ≥ 0. In addition, particles can replicate and die with state-dependent rates. This leads to the following re-parameterization of the generic model defined above: The population densities are denoted as n + and n − , and their dynamics is described by: The front velocity is obtained by substituting equation (13) into equation (11).
Compared to the simple reaction-diffusion-advection equation, persistent random walk model has more parameter. Instead of advection velocity and diffusion constant, there are two velocities and two switching rates. Furthermore, there are four replication rates instead of one. Below we describe how each of these parameters influences front propagation and compare the results to the prediction of the simpler reaction-diffusion-advection model equation (1) (to which we refer to as the diffusion approximation).
Most of the discussion below is based on the exact solution. To check its validity, we performed a few numerical simulations (see Appendix A 6 for details). Numerical results are presented in figure 3 to 5, which show excellent agreement between the theory and simulations.
Given the large number of parameters it is important to develop some intuition first. To this end, we begin by considering a simple model with a single replication rate r ++ > 0, and all other replication and death rates equal to zero. This simplified case is used to illustrate the effects of persistence on front propagation and compare the exact solution and the diffusion approximation. Then, we return to the most general case and determine how various elements of the replication matrix influence front propagation. Finally, we analyze the effect of death rates.

A. Comparison between the exact solution and the diffusion limit
With only r ++ not equal to zero, the expression for the front velocity simplifies to When r ++ = 0, the velocity is given by the time-averaged velocity of the random walkerv. This average velocity also equals to the advection velocity U in the diffusion approximation.  (14). Particles interconvert between two states: a "+ state" that moves to the right at speed v+ and a "− state" that moves to the left at speed v−. The transition frequencies are denoted by f+ and f−. Note that a minimal replication rate is required for the front to propagate rightward whenv < 0 and that V = v+ when r++ ≥ f+.
To avoid confusion and simplify the notation, we refer to both of these velocities asv. The detailed derivation of the diffusion approximation is provided in Appendix A 4.
The major difference between the exact solution and the diffusion approximation is their dependence on r ++ . In the diffusion approximation, the front velocity grows indefinitely with r ++ , but the exact solution correctly recapitulates the physical constraint that V ≤ v + . We illustrate this difference in figure 3 where we show V (r ++ ) for three parameter combinations that correspond tov > 0,v = 0, andv < 0. Diffusion approximation works best forv near zero and small replication rates. Deviations from this sweet spot lead to drastic differences, and the diffusion approximation can both overestimate (figure 3a) and underestimate (figure 3c) the velocity of the front.
In contrast to replication, the dependence of V on v + and v − is largely determined by dimensional analysis and Galilean invariance (see Appendix A 3). In particular, V must be a linear combination of v + and v − of the following form:  (2)) and our exact solution with solid lines (equation (15)). The simulation results are shown with circles. The effective parameters in the diffusion approximations are derived in Appendix A 4.
Both the exact solution and the diffusion approximation comply with equation (16) and therefore should show similar dependence on the microscopic velocities. figure 4 confirms this expectation. When the ratio of v+ v− is held fixed (figure 4a), the expansion velocity is linear in v + + v − for both equation (2) and (15) although the slopes of the lines do differ. When the sum of the velocities is held fixed (figure 4b), the dependence on the velocity ratio is also linear, but now the diffusion approximation and the exact solution differ only in the intercept, again consistent with equation (16).
The effects of switching frequencies are somewhat similar to those of r ++ because V depends only on f + /r ++ and f − /r ++ . Both the exact solution and diffusion approximation predict that the front velocity decreases with increasing the total switching frequency (figure 5a) and the relative magnitude of f + (figure 5b). The values of the velocity can, however, be very different. We find the largest deviations between the diffusion approximation and the exact solution for small switching frequencies. In this regime, the diffusion approximation predicts velocities above v + just as for large replication rates (compare to figure 3).
Finally, we show two qualitatively different front shapes in figure 6. For r ++ < f + , the front has sigmoidal shape with exponential decay at the leading edge. Qualitatively, this is the same as in the diffusion approximation. For r ++ > f + , the leading edge of the front is infinitely sharp, i.e. it is effectively a step function followed by gradual saturation towards the carrying capacity. The origin of this singularity is that the front velocity saturates at v + and the shape of the leading edge does not relax, but preserves any discontinuities present in the initial condition (we start simulations with a nonzero density within an interval).

B. Contribution of different replication modes to front propagation
In the preceding section, we considered replication only in the plus state. Now, we systematically compare how different modes of replication contribute to the onset of propagation For small replication rates, we find that V αβ ≈ U + C α,β √ r α,β ; see Appendix A 5 a. The constants C α,β coincide for the V ++ and V −+ and for V +− and V −− . This clustering into two pairs is evident both in figure 7a and in figure 7b.
In the opposite limit of high replication rates, the front velocity reaches the physically maximal value of v + (see Appendix A 5 a). While V +− , V −− , and V −+ approach v + gradually as r αβ → ∞, V ++ becomes equal to v + when r ++ = f + . In the former cases, right moving particles need to pass through left-moving state during replication, so V < v + because some time is spent in the left-moving state. In the latter case, however, some particles always stay in the right-moving state when the production rate r ++ exceeds the loss rate f + , so they create a front moving with V = v + .
When U < 0, front velocities have different sign for small and large replication rates. The condition V α,β = 0 marks the onset of propagation strictly to the right, which is important in many applications 21 . Below we provide critical replication rates r * αβ necessary for a right-moving front.
• only r ++ = 0 • only r +− = 0 • only r −− = 0 • only r −+ = 0 In all cases, parameters that favor front propagation in the positive direction lower the critical replication rate.
Finally, we determine how the ordering of the relative magnitudes of V α,β as well as r * αβ depend on model parameters. In Appendix A 5 b, we analyse all six pairwise comparisons of the four velocity equations. We first find that the ordering of front velocities depends the dimensionless parameter µ = v+ v− , which defines regimes A, B, and C with unique ordering and crossing of velocities. The ordering of critical replication rate additionally depends on the value of ν = f− f+ , which leads to additional possibilities labeled by α, β, γ, δ, . In total, there are eight possible regimes: Aα, Bβ, Bγ, Bδ, B , Cγ, Cδ, and C (see figure 10 and Table I in Appendix A 5 b). These results are summarized visually by a phase diagram in figure 7b.

C. Particle death leads to a discontinuous jump in the front velocity
In this section, we show that nonzero death rates could lead to a qualitatively new behavior: the system can transition abruptly from extinction to a traveling wave moving at a finite velocity ( figure 8). This transition occurs when the net growth rate is positive in one, but negative in the other state.
The transition between extinction and growth does not depend on particle velocities. Indeed, let us integrate the linearized equation (5) over space: The boundary terms disappear because the concentration is zero at ±∞, and we obtain the non-spatial version of the model: where C α (t) = ∞ −∞ n α (t, x)dx.
Extinction (C α = 0) occurs only when all the eigenvalues of Λ α,β have negative real part.
For the two-species system, the transition from extinction to growth occurs when one of the eigenvalues crosses zero. Thus, generically, we expect the transition when the trace is negative Λ 22 + Λ 11 < 0 and the determinant is zero: We illustrate the transition to growth further using the following simple, yet generic, scenario with r ++ > 0, d − > 0, and r +− = r −+ = d + = r −− = 0. The minimal replication rate r gap + is given by Note that r gap + is different from the critical replication rates equation (17)- (20) for which the front velocity is zero. At the transition, the velocity of the front is given by Following Ref. 32 , we term this minimal velocity as the gap velocity since it characterizes the gap in the possible values that front velocity can take. The diffusion limit predicts both r gap + and V gap correctly; see equations (A36) and (A37) .

V. ORIGINS OF GAP VELOCITY IN REACTION-TRANSPORT SYSTEMS
Onset of growth with a non-zero velocity has been recently reported in a growing network of microtubules 32 . Based on their model, the authors of Ref. 32 argued that the velocity gap is a unique feature of persistent motion with replication. Here, we re-examine this claim and identify the key differences between populations of persistent random walks and persistentlygrowing polymers. We first review the polymer model and compare its behavior to diffusion approximation and replicating particles that do not die. We then revisit this analysis by including death rates to account for the possibility of complete polymer depolymerization.

A. Propagating fronts of replicating polymers undergoing persistent growth dynamics
We begin by briefly summarizing the key results of the model for replicating polymers undergoing persistent growth 32 . In this model, polymers have a static that does not move and a dynamic end at which the polymer can polymerize with velocity v + or depolymerize with velocity v − . It is assumed that the dynamic end is to the right of the static end due to the mechanical interactions between the polymers in the network. This implies that, if we define the position of the static end as x s and the polymer length as l, then the position of the dynamic end is x d = x s + l. We denote the polymerizing or growing state of the dynamic end as + and the depolymerizing or shrinking sate as −. The transition rate from + to − sate is f + . The rate of the reverse transition is f − . Replication in this model is the nucleation of a new zero-length polymer in the growing sate. We denote this replication rate by Q (t,x) .
The mathematical formulation of the above-described polymer model reads where ρ + (t, x, l) and ρ − (t, x, l) are the densities of polymers in growing and shrinking states respectively at time t with their static ends located at x and length equal to l. We use δ(l) to denote the Dirac delta function.
The basic model of polymer replication considered in Ref. 32 assumes that each growing end can branch. This makes Q a function of only the density of the growing ends, C + (t, x + ), which can be computed from ρ + (x s , l) by taking an integral over x s and l while keeping x s + l = x; see Ref. 32 . Upon assuming simple logistic-like saturation at high C + , we obtain the following equation for Q where K controls the saturation density, and R + is the replication rate.
The calculations in Ref. 32 show that this system can either become extinct or expand as a traveling wave. Note that the wave travels in x-space with the distribution of polymer lengths being constant at each spatial location in the co-moving reference frame. The velocity of the front is given by and the critical replication rate is At the transition from extinction to expansion, the front velocity jumps from zero to the gap velocity

B. Comparison between propagation onset in persistent random walks and polymers
The theory of persistent polymers predicts a minimal spreading velocity given by equation (30). Is this a qualitatively new feature of persistent polymers or can this transition be captured by persistent random walk model or even the diffusion approximation?
To answer this question, we first compared the predictions of the polymer model and the persistent random walk model with identical parameters. Typical dependence of the velocity on the replication rate is shown in figure 9a. As suggested in Ref. 32 , there is indeed a qualitative difference between the two models. Propagation begins with zero velocity in the random walk model, but with a nonzero gap velocity in the polymer model.
Given that we observed gap velocity for the random walk model with death, we hypothesized that the discontinuous jump in the polymer model is related to death-like events. Indeed, when polymer shrinks to zero length, it disappears from the population, i.e. it effectively dies. The effective death rate can be estimated as the inverse lifetime of the polymer, which has been computed in Ref. 32 . Using their results, we then investigated the persistent random walk model with d + = 0 (growing polymers never die) and The results of this comparison are shown in figure 9b. Now, both the polymer and the random walk model have indistinguishable values of the critical replication rate and the gap velocity. This agreement holds for all parameter values. Indeed, one can confirm this by substituting the death rate from equation (31) into equation (25) for the gap velocity in the random walk model. The result is identical to equation (30) for the gap velocity in the polymer model. Thus, the sudden onset of growth in the polymer model can be traced to the disappearance of polymers that shrink to zero length. The persistent nature of polymer growth is less important because the diffusion approximation predicts identical values of the critical nucleation rate and the gap velocity.
figure 9b also shows that the agreement between the random walk and polymer models is only qualitative. Both reproduce the discontinuous onset of growth and saturation of the front velocity at v + , but the intermediate values of the front velocity do differ. This difference is however very small for many biologically relevant parameter values. Therefore, one might be able to use a much simpler and substantially more computationally efficient random walk model to approximate the dynamics of polymer networks.

VI. DISCUSSION
Front propagation in replicating and persistently moving particles is a widespread phenomenon that has many applications 1,15 . Such reaction-transport systems can be approxi-  (15)), and the self-replicating polymer model (dashed and dotted lines, equation (28)). Only the polymer model has a jump in front velocity when the replication rate is increased. We used the parameters of the polymer model to parameterize the other models: When death is incorporated using equation (31), the polymer and the random walk models show qualitatively the same behavior.
mately mapped onto well-understood reaction-diffusion systems, but the diffusion approximation can be rather inaccurate or even produce unphysical results. Furthermore, the diffusion approximation is insensitive to certain parameters of the persistent random walks that do affect the behavior of the system.
In this paper, we formulated a general model of persistent random walks with an arbitrary number of motility and replication states, and provided a simple method to compute the expansion velocity. The presented framework is very flexible and can account for complex transitions between states. For example, a time delay could be modeled by requiring that any transition occurs through one or more intermediate states.
The general framework was then applied to a two-state system. We obtained an exact solution for the most general model and explored the dependence of the front velocity on all model parameters. We found that the diffusion approximation is accurate only for a limited range of parameter values. In particular, the diffusion approximation fails for high replication or low switching rates ( figure 3 and 5).
Furthermore, the diffusion approximation does not fully capture the dependence of the front velocity on replication rates in different states. We analyzed the effect of four possible replication mechanisms: two choices for the state in which the particle replicates and two choices of the offspring state. These replication modes make unequal contributions to the front velocity, and the efficacy of each replication mode depends on other parameters. We completely characterized these dependencies and identified eight regions in the parameter space that exhibit different dependence of front velocities on replication rates; see figure 7.
We also analyzed the effects of particle death, which have been rarely considered in the literature. The death rates result in somewhat unusual behavior: There is a minimal speed at which the front can propagate. This minimal speed occurs right when the replication can counterbalance death, and the population starts to grow. Such transitions from no growth to rapidly moving expansion fronts have been observed in microtubule networks 32 and we sought to provide an intuitive explanation for their origin. We found that one can identify an effective death term in the persistent polymer model, which fully accounts for the critical replication rate and velocity jump. In fact, the persistent particle model approximates the persistent polymer model rather well across the whole range of replication rates (figure 9b).
Overall, our work extends earlier results on persistent random walks 1-6 to a much more general class of models and provides new insights into the role of different replication modes and death rates. These results might find applications across a number of fields where one encounters front propagation. Our goal is to compute the velocity of traveling fronts described by equation (5), which we repeat below

ACKNOWLEDGMENTS
First, we perform the Fourier transform in the space domain, x → k: Then, we perform the Laplace transform in the time domain t → s: The last equation can be recast in the following form β [(s + iv α k)δ αβ − Λ αβ ] n β (s, k) = n α (t = 0, k), which immediately suggests an equivalent matrix formulation: where A αβ = (s + iv α k)δ αβ − Λ αβ . The solution of the matrix equation reads We now perform the inverse Laplace transform, s → t: The residues arise when detA = 0. The degeneracy of A implies that there exist h = 0 such that Ah = 0. Therefore, sh α = (−iv α kδ αβ + Λ αβ )h β , i.e. s is an eigenvalue of B αβ = Λ αβ − iv α kδ αβ . Generically, B has N eigenvalues s (l) , which define distinct branches of the dispersion relations with corresponding velocities V (l) . We expect that the branch with the largest front velocity describes the long-time behavior of our system. Next, we perform the inverse Fourier transform, k → x: [n α (s, k)] . (A8) In the above derivation, we assumed the generic situation of distinct eigenvalues of B, in which case all poles are first order. In the middle line, we have transformed the spatial variable to the reference frame comoving with the front: ξ = x − V t, where V is the front velocity we wish to determine.

b. Asymptotic evaluation of the inverse Fourier transform
The integrals in equation (A8) can be evaluated in the long time limit using the steepest descent or saddlepoint method; see Ref.
2 . The controlling factor is the exponential term, so we are looking for the saddlepoint of s(k) + ikV . We also impose time invariance of the front by requiring that the real part of s(k)t + ikV equals zero; otherwise, the particle density in the comoving frame either diverges or vanishes. Thus, we obtain where k f is the value of k at the saddle point. It is easy to see that k f is purely imaginary, so we let k f = iκ f . Then, the system of equations that we need to solve takes the following form where all variables are real. The solution of these equations defines the front velocity for each branch l. The actual velocity is obtained by making the largest of l values.
In the following, we will first solve for κ f by eliminating V from the system of equations, which yields the following equation for κ f : 2. Front velocity for the two-state persistent random walk We start from the compact formulation of the problem and quickly repeat the analysis outlined in the preceding section: Fourier transform x → k for the spatial coordinate and the Laplace transform t → s for the temporal coordinate yields: sn 1 − n 1 (t = 0, k) =−ikv 1 n 1 + Λ 11 n 1 + Λ 12 n 2 sn 2 − n 2 (t = 0, k) =−ikv 2 n 2 + Λ 21 n 1 + Λ 22 n 2 .
We recast this result in the matrix form where We obtain s(k) by setting the determinant of A to zero Instead of directly solving the above equation and substituting it into equation (A11), we differentiate equation (A16) first, then solve for ds/dκ, and finally substitute the result into equation (A11). This procedure yields additional equation We now need to solve equation (A17) and equation (A16) simultaneously. This is can be done by constructing a linear combination of these equations to eliminate s 2 and then solving a linear equation for s in terms of κ. The result can be substituted in either equation and solved to obtain κ. The results read The front velocity is then obtained from equation (A10). By examining the cases for Λ 11 + Λ 22 < 0 and Λ 11 + Λ 22 > 0 separately, we find that, for either case, the front velocity that corresponds to the right side of the population is The persistent random walk model must be invariant under Galilean transformation. We can use this invariance to determine how the front velocity depends on model parameters. Indeed, consider shifting the reference frame from the lab frame to the frame moving with velocity u. In other words, v i → v i − u and V → V − u. Then, we expect that Differentiating this with respect to u, we find Solving this equation for the two-state system gives where w is a function that depends only on the switching, replication, and death rates. Thus, Galilean invariance completely determines how the front velocity depends on the velocities of the states. Furthermore, equation (A25) is also very useful for simulations because one can, for example, move into a reference frame in which the velocity in the two states are equal in magnitude, but opposite in orientation. Such a choice substantially improves the accuracy of fine difference methods; for further details see Appendix A 6 on numerical simulations.

Diffusion approximation
Here, we show that the classical reaction-advection-diffusion equation (1) and its front velocity (2) are a limiting case of the persistent random walk model. We consider a twostate system with right and left moving states as in equation (14). The advection velocity, or the time-weighted average velocity in the lab reference frame, is given bȳ To simplify the derivation, let us change into a reference frame that moves with velocity u = (v + − v − )/2. Then, the velocities of left and right moving states have equal magnitude v = (v + + v − )/2, and equation (14) takes the following form which via the relations (13) is equivalent to From the first equation in (A28), we obtain Upon substituting this in the second equation in (A28) and rearranging the terms, we find that where ω = −Λ 11 − Λ 22 . We now assume slow temporal evolution t ∼ x 2 and neglect ∂ 2 ∂t 2 term, then equation (A30) reduces to equation (1) with where U u is U in the reference frame moving with velocity u, and Note that the two terms in the approximate expression for g are the time averaged replication and death rates.
We obtain the final result by returning to the lab reference frame, which gives and the front velocity of the right side of the population The critical replication rate corresponding to the transition from extinction to expansion occurs when g, defined by equation (A33), crosses zero. This is given by the condition which is identical to equation (23) for persistent random walks.
The corresponding gap velocity is given by equation (A35) at g = 0. Hence, we have which matches the corresponding expression for persistent random walks.

Front propagation with a single mode of replication a. Front velocity solutions V αβ
Here, we provide the front velocity solutions when only one replication rate r αβ is nonzero. We denote V (r ++ ≥ 0, r +− = r −+ = r −− = 0) as V ++ and similarly for other replication modes.
The behavior for small r ++ is as follows is the velocity at zero growth and death rates as given by equation (12).
For r ++ ≥ f + , the profile becomes infinitely sharp, and the velocity equals v + . As r ++ approaches f + from below, the velocity approaches v + quadratically: • r +− only For small r +− , For large r +− , the velocity approaches v + as For small r −− , • r −+ only For small r −+ , For r −+ → ∞, the velocity approaches v + as Ordering of the front velocities V αβ Here we examine the relative magnitudes of four velocities V αβ arising from a single mode of replication. From the front velocity equations (A38), (A41), (A44), (A47), we perform a pairwise analysis and ask the following: • Is one velocity V αβ always greater than the other? Or, as the replication rate is increased, does one velocity surpass the other with increasing r αβ ?
• Given that one velocity surpasses the other, what is the condition for one critical replication r * αβ to be larger than that of the other?
The answer to the first question depends on the dimensionless parameter ν = f− f+ , while that to the second question additionally depends on µ = v+ v− . We denote the non-dimensionalized replication rates as ρ = r αβ f+ to directly compare different modes of replication.

Numerical simulations
We simulated the continuum equations using an explicit first order finite difference method.
To avoid numerical diffusion, the exact solution of the advection equation was used for the numerical updates. Time and space discretizations were chosen to ensure v α ∆t ∆x was a rational number, and concentrations of each type were periodically shifted by the appropriate integer number of lattice points. Simulations were initialized as a step function in the middle of the of the lattice. We imposed absorbing boundary concentrations at the boundaries, and ensured that simulations ended before they reached either boundary. The velocity was obtained from fitting the front position, defined as a point with concentration closest to 10 −3 , to a linear function of time.