An optimal control approach to nonlinear wave speed selection in reaction-diffusion equations

An optimal con trol approac h to nonlinear w a v e sp eed selection in reaction–diffusion equations Reb ecca M. Crossley 1 , Carles F alcó 1 , and Ruth E. Baker 1 1 Mathematic al Institute, University of Oxfor d, Oxfor d, Unite d Kingdom OX2 6GG T rav elling w av e solutions of reaction–diffusion equations are widely used to mo del the spatial spread of populations and other phenomena in biology and physics. In this article, w e rein terpret the classical v ariational principle approac h through an optimal control form ulation, in order to obtain a lo wer b ound on the inv asion sp eed of trav elling w av e solutions in systems of nonlinear partial differen tial equations. W e b egin b y analysing single–sp ecies mo dels, where the evolution of the density is go v erned b y a scalar equation with a density-dependent diffusion term and a nonlinear reaction term. W e show that for an y admissible test function, maximising with resp ect to the parameter of interest yields a b ound on the trav elling wa v e sp eed. W e apply this framework to several examples, including the p orous–Fisher equation, and examine when nonlinear selection mechanisms dominate o ver the classical linear marginal stability criterion. Extending this approach, we then consider m ulti-sp ecies systems of reaction–diffusion equations and, reframed as Pon try agin-type optimalit y systems, we deriv e analogous b ounds on the trav elling wa v e sp eed using a v ariational framework under weak coupling. Finally , w e employ numerical sim ulations to confirm the accuracy of the predicted wa ve sp eeds across a range of illustrative examples. Key words. trav elling wa v es, nonlinearity , optimal control, v ariational principle, reaction–diffusion MSC co des. 35A23, 35C07, 49J40, 49K20, 92-10 1 In tro duction Reaction–diffusion equations ha ve long served as a p o w erful and widely-used mathematical framework for mo delling the collectiv e mo vemen t of p opulations of individuals. These mo dels typically tak e the form of partial differential equations (PDEs) with density-dependent reaction (growth) and diffusion terms, enabling a detailed description of b oth the spatial and temp oral ev olution of a p opulation across a v ariet y of biological and ecological contexts [1]. A hallmark feature of such models is the emergence of trav elling wa v e solutions that maintain a fixed profile while propagating at a constan t sp eed. These trav elling w av es hav e b ecome central to our theoretical understanding of diverse biological pro cesses, from w ound healing to tumour in v asion, and sp ecies migration [ 2 – 6 ]. In many of these settings, the wa v e connects an unstable steady state ( e.g. , uncolonised space) to a stable steady state ( e.g. , an established p opulation), even in the absence of 1 classical linear instabilities. This is unlike systems such as crystal growth, nerve impulses, or fluid flo w, where instabilities in the linear system are key driv ers of wa v e propagation. One of the seminal mo dels exhibiting trav elling wa v e solutions is the Fisher–KPP mo del [ 1 , 7 , 8 ]. Along with its many extensions, the Fisher–KPP mo del clearly illustrates ho w changes in the mo del parameters shap e the sp eed and form of the inv ading wa v e [ 1 , 7 – 9 ]. Considerable analytical effort has fo cused on understanding the influence of the v arious nonlinearities in single–sp ecies reaction– diffusion equations, with recen t w ork yielding increasingly refined existence and predicted wa ve sp eed results [4, 10–19]. Sev eral complementary p ersp ectiv es ha ve emerged to explain ho w trav elling w av e sp eeds are selected, eac h with its o wn strengths and limitations. A foundational viewpoint uses the linear marginal stabilit y criterion, whic h identifies the trav elling wa v e with the steep est p ermissible exp onen tial decay in to the unstable state [ 20 ]. While this criterion provides an elegant and computable prediction, it relies on linearisation around the leading edge and therefore cannot accoun t for cases where nonlinear in teractions—either in the bulk or in the front—influence the selected sp eed. T o address this limitation, Benguria and Depassier introduced a v ariational metho d for wa ve sp eed selection, initially for nonlinear reaction terms [ 21 , 22 ] and later extended to incorp orate nonlinear diffusion [ 23 ]. Their approac h generates rigorous upp er and low er b ounds on the admissible w av e sp eeds and, crucially , pro vides a selection criterion identifying when the minimal admissible sp eed is indeed the realised one [ 24 ]. Hadeler and Rothe [ 15 ] clarified that, under suitable conditions, the asymptotic propagation sp eed can coincide with the linearly predicted v alue, ev en when the v ariational upp er and lo wer b ounds do not match. These results highligh t cases where the linear and nonlinear sp eed selection criteria actually coincide. Though broadly undev elop ed thereafter, more recent studies, suc h as that of Stok es et al. [ 25 ], hav e since extended these ideas by systematically analysing how densit y-dep enden t diffusion shap es b oth tra velling wa v e existence and sp eed, thereb y rev ealing parameter regimes where nonlinear diffusion either reinforces or inv alidates predictions from earlier linear or v ariational theories. T aken together, these developmen ts establish a ric h set of to ols for analysing trav elling w av es in single–sp ecies systems, y et they also exp ose a notable gap: most classical sp eed-selection principles, including b oth the marginal stability criterion and the Benguria–Depassier v ariational framework, are form ulated for scalar equations. How ever, man y biological inv asion pro cesses fundamen tally inv olve in teractions b et w een m ultiple p opulations. Extending the theory to determine trav elling wa ve sp eed selection in these multi-species systems is therefore essential for connecting mechanistic mo dels of collectiv e migration with the analytical insights that hav e prov ed so v aluable in the single–sp ecies setting. W e analyse a range of mathematical mo dels with differen t combinations of densit y-dep enden t reaction and diffusion terms. Analytical tec hniques are used where feasible, supp orted by n umerical sim ulations implemen ted using a finite-difference sc heme describ ed in [ 9 ], to illustrate the existence and selection of trav elling wa v es across differen t parameter regimes. The main con tributions of this pap er are threefold. First, we rein terpret the classical v ariational principle for trav elling w av e sp eed selection as an optimal control problem, clarifying its structure and limitations. Second, we extend this formulation to multi-species reaction–diffusion systems, where classical v ariational principles are 2 no longer directly applicable. Third, we derive a tractable v ariational appro ximation of the resulting optimalit y system that yields sharp low er b ounds on trav elling wa v e sp eeds in a range of biologically motiv ated mo dels, including systems exhibiting nonlinear sp eed selection. The remainder of the pap er is organised as follows. First we presen t the optimal control approac h to nonlinear wa ve sp eed selection for general multispecies mo dels (Sec. 2). Next, in Sec. 3 we revisit the classical v ariational principle for a single–sp ecies and pro vide a sufficien t criterion for nonlinear w av e sp eed selection in terms of the functional forms of the diffusion and reaction terms. This scalar setting provides a b enc hmark against which we can assess how sp eed-selection mechanisms change when additional in teractions are incorp orated. W e then apply the optimal control framew ork to a n umber of examples in the literature including for mo dels inv olving t wo in teracting sp ecies (Sec. 4). W e show that this approac h successfully predicts a minimum trav elling wa ve sp eed in a n umber of examples, including cell inv asion and differentiation. W e conclude in Sec. 5 with a brief discussion of these results and suggest some p otential av enues for future research. 2 An optimal control approac h for m ultisp ecies models W e consider reaction–diffusion systems for n ≥ 1 sp ecies ρ = ( ρ 1 , ˜ ρ ) ⊤ , with ˜ ρ = ( ρ 2 , . . . , ρ n ) ⊤ , of the form ∂ t ρ 1 = ∂ x ( D ( ρ ) ∂ x ρ 1 ) + f ( ρ ), (1) ∂ t ˜ ρ = g ( ρ ), (2) where ρ i = ρ i ( x , t ) with x ∈ R , t ≥ 0 , and D , f : R n → R , g : R n → R n − 1 . This type of mo del considers a single motile p opulation, ρ 1 , with the rest of sp ecies ev olving via reaction terms—these are often used, for instance, in systems describing cell inv asion into the extracellular matrix [ 5 , 9 , 11 ]. W e fo cus on systems admitting constant sp eed ( c ≥ 0 ), constant profile trav elling wa v e solutions ρ ( x , t ) = u ( z ) , z = x − ct , where u = ( u 1 , ˜ u ) ⊤ and ˜ u = ( u 2 , . . . , u n ) ⊤ . W e assume that u 1 is non-increasing and, after nondimensionalisation,    u 1 ( z ) → 1 as z → −∞ , u 1 ( z ) → 0 as z → + ∞ . The trav elling wa ve solution also fixes b oundary conditions for ˜ u for z → ±∞ . In trav elling wa ve co ordinates, and defining v = − d u 1 / d z ≥ 0 , the system b ecomes d d z ( D ( u ) v ) = − cv + f ( u ), c d ˜ u d z = − g ( u ). Using the monotonicity of u 1 , we can rewrite the system in terms of u 1 , with d / d z = − v ( u 1 ) d / d u 1 , 3 yielding v ( u 1 ) d d u 1 ( D ( u ) v ( u 1 )) = cv ( u 1 ) − f ( u ), (3) cv ( u 1 ) d ˜ u d u 1 = g ( u ). (4) T o obtain a v ariational characterisation of the w av e sp eed c , w e introduce a non-negative, non- increasing test function φ ( u 1 ) , with ϕ ( u 1 ) = − φ ′ ( u 1 ) ≥ 0 and φ (1) = 0 . W e note that for systems ( n > 1 ), the remaining equations (Eq. (4) ) act as dynamical constrain ts. Multiplying Eq. (3) b y D ( u ) φ ( u 1 ) and applying in tegration b y parts gives I c [ v ; ˜ u ] : = Z 1 0  − 1 2 ϕ ( u 1 ) D ( u ) 2 v ( u 1 ) 2 + cv ( u 1 ) D ( u ) φ ( u 1 ) − D ( u ) f ( u ) φ ( u 1 )  d u 1 = 0. The quadratic dep endence of I c on v suggests studying the asso ciated v ariational problem, with the aim of deriving b ounds that may pro vide direct estimates of the wa v e sp eed, c , via 0 = I c [ v ; ˜ u ] ≤ I c [ v ∗ ; ˜ u ∗ ] , for an optimal pair ( v ∗ , ˜ u ∗ ) . In systems, how ev er, v app ears implicitly through Eq. (4) . W e therefore consider the optimisation problem v ∗ = arg max v I c [ v ; ˜ u ] , sub ject to cv ( u 1 ) d ˜ u d u 1 = g ( u ). (5) This is a standard optimal control problem, which we solv e using Pon try agin’s maxim um princi- ple [26]. Introducing the Hamiltonian defined by H [ u 1 , v ; ˜ u , ω ] = − 1 2 ϕ ( u 1 ) D ( u ) 2 v ( u 1 ) 2 + c v ( u 1 ) D ( u ) φ ( u 1 ) − D ( u ) f ( u ) φ ( u 1 ) + ω ⊤ · g ( u ) cv ( u 1 ) , where ω = ( ω 2 , . . . , ω n ) ⊤ are co-state v ariables. Pon tryagin’s principle implies that maximisers v ∗ are p oin t wise maximisers of H . When the Hamiltonian is concav e down in v , these satisfy ∂ H ∂ v = − ϕ ( u 1 ) D ( u ) v ( u 1 ) + c D ( u ) φ ( u 1 ) − ω ⊤ · g ( u ) cv ( u 1 ) 2 = 0 , or equiv alently − ϕ ( u 1 ) D ( u ) 2 v ( u 1 ) 3 + c D ( u ) φ ( u 1 ) v ( u 1 ) 2 − ω ⊤ · g ( u ) c = 0 . (6) Eq. (6) implicitly determines the optimiser v ∗ as the positive real ro ot of the cubic equation. 4 F urthermore, the co-state v ariables evolv e according to d ω d u 1 = −∇ ˜ u H =  ϕ ( u 1 ) D ( u ) v ( u 1 ) 2 − cv ( u 1 ) φ ( u 1 )  ∇ ˜ u D ( u ) + ∇ ˜ u ( D ( u ) f ( u )) φ ( u 1 ) − ω ⊤ ∇ ˜ u g ( u ) cv ( u 1 ) . No transv ersality conditions are imp osed on ω , since trav elling w av e solutions prescrib e b oundary conditions for the state v ariables ˜ u at u 1 = 0 and u 1 = 1 . In this next section, w e analyse this problem for the sp ecific case of one-sp ecies mo dels. In doing so, w e unv eil a relationship to the v ariational principle for one-sp ecies models, previously analysed in [ 21 , 22 , 24 , 25 ], b efore applying the optimal con trol approach to a num b er of new examples that illustrate why the optimisation yields useful b ounds on the tra velling w av e sp eed c . 3 One-sp ecies mo dels: reco vering the v ariational principle In cases where there is only one species, ρ = ( ρ 1 ) , the dynamical constrain t imp osed by Eq. (4) v anishes. Here, we require that f ( ρ 1 ) has tw o equilibrium p oints, we choose these without loss of generalit y at ρ 1 = 0 , and ρ 1 = 1 , i.e. , f (0) = f (1) = 0 . In this case, maximising the functional I c relies on standard v ariational argumen ts. In particular, omitting indices for simplicity , we hav e I c [ v ] := Z 1 0  − 1 2 ϕ ( u ) D ( u ) 2 v ( u ) 2 + c v ( u ) D ( u ) φ ( u ) − D ( u ) f ( u ) φ ( u )  d u = 0. (7) In terpreting the integrand in Eq. (7) as a quadratic p olynomial in v reveals a natural inequality: − 1 2 ϕ ( u ) D ( u ) 2 v ( u ) 2 + c v ( u ) D ( u ) φ ( u ) ≤ c 2 φ ( u ) 2 2 ϕ ( u ) , (8) that transforms the in tegral identit y into a b ound that holds for all admissible test functions. Equalit y in Eq. (8) is only obtained when v ( u ) = v ∗ ( u ) := cφ ( u ) D ( u ) ϕ ( u ) . (9) Inserting the b ound in Eq. (8) in to Eq. (7) yields a family of low er b ounds on the wa ve sp eed, indexed b y the choice of smo oth test function φ ( u ) with φ ′ ( u ) = − ϕ ( u ) : 1 2 c 2 ≥ R 1 0 D ( u ) f ( u ) φ ( u ) d u R 1 0 φ ( u ) 2 /ϕ ( u ) d u , (10) T o make the b ound explicit, w e now select a conv enient one-parameter family of test functions, φ β ( u ) = ((1 − u ) /u ) β , for β ∈ [0, 2) . This choice, previously used in [ 24 , 25 ]; satisfies the condition that the integral remains finite and b oundary terms v anish in the limiting sense. As such, it balances analytical tractabilit y with the ability to capture the asymptotic b eha viour near the steady states. 5 Using this test function, we find 1 2 c 2 ≥ sup β ∈ [0,2) F ( β ), with F ( β ) = β R 1 0 D ( u ) f ( u ) u − β (1 − u ) β d u B (2 − β , 2 + β ) , (11) where B is the Beta function. T o confirm that the inequality in Eq. (11) defines a genuine v ariational principle, we must also sho w that there exists a test function, φ = ˆ φ (and corresp ondingly , ϕ = ˆ ϕ ) for which the b ound is attained. This requires analysing the asymptotic b eha viour of the trav elling wa v e profile near u = 0 and u = 1 and verifying that the corresp onding integrals con verge. Equality in Eq. (10) is attained (as p er Eq. (9)) when v ( u ) = c ˆ φ ( u ) D ( u ) ˆ ϕ ( u ) = − c ˆ φ ( u ) D ( u ) ˆ φ ′ ( u ) , (12) where ′ = d / d u . Eq. (12) can b e integrated to obtain ˆ φ ( u ) = exp ( − c Z u ¯ u d q D ( q ) v ( q ) ) , (13) for some ¯ u ∈ (0, 1). Clearly ˆ φ > 0 and is monotonically decreasing. W e no w make the c hange of v ariables using p ( u ) = D ( u ) v ( u ). Eq. (3) b ecomes p ( u ) d d u p ( u ) − cp ( u ) + D ( u ) f ( u ) = 0. (14) In order to ensure conv ergence of ˆ φ ( u ) across all of u ∈ [0, 1] , we must then analyse the asymptotic b eha viour of Eq. (14) as u → 0 + and u → 1 − [27]. Begin by considering the linear analysis around the endp oin t u = 0. Near u = 0 , we can write f ( u ) ∼ f ′ (0) u and D ( u ) ∼ D (0) (if D (0) = 0 or f ′ (0) = 0 , then the marginal stability analysis predicts a zero wa ve sp eed), so that b y emplo ying the ansatz p ( u ) = αu , where α is the larger ro ot of m 2 − cm + D (0) f ′ (0) = 0 , we find that c ≥ 2 p D (0) f ′ (0) = c L , in order for α to b e real, and hence av oid oscillatory solutions that predict negative densities. Hereon in, c L is the trav elling w av e speed predicted by linear stabilit y analysis. Solving Eq. (13) for p ( u ) = αu and ¯ u = u ( z → −∞ ) = 1 giv es ˆ φ = u − c/α , which ensures that all in tegrals in Eq. (10) con verge as long as c/α < 2. No w consider the linear analysis around the endpoint u = 1 . Define the new v ariable, s , as s = 1 − u . Using T aylor series expansions, f ( u ) ∼ f ′ (1)( u − 1) = −| f ′ (1) | s and D ( u ) ∼ D (1). Setting p ( u ) = r (1 − u ) = r s and employing this ansatz in Eq. (14) then, following the same metho ds as b efore, w e find that for ¯ u = 1 , ˆ φ = (1 − u ) c/r with r as the p ositiv e ro ot of the auxiliary equation 6 r 2 + cr + f ′ (1) D (1) = 0 , giving r = − c ± p c 2 − 4 D (1) f ′ (1) 2 . As such, the integrals in the functional (11) m ust conv erge at u = 1 b ecause c/r > 0 and do es not restrict the sp eed. This is sufficien t to sho w conv ergence of all of the integrals across the entire domain [0, 1] for an attainable test function ˆ φ , so that Eq. (11) do es in fact define a v ariational principle. 3.1 A criterion for nonlinear w av e sp eed selection. The classical v ariational framew ork pro vides low er b ounds on the admissible trav elling wa ve sp eed but do es not, in general, indicate whether the minimal sp eed is selected b y linear leading-edge dynamics or b y nonlinear effects in the in terior of the wa v e. In this section, w e extract from the v ariational form ulation a simple diagnostic criterion that predicts when nonlinear (“pushed”) wa ve sp eed selection o ccurs, i.e. , when c > c L . The criterion is based on the b eha viour of the v ariational functional F ( β ) (Eq. (10) ) near the b oundary β = 2 , which corresp onds to the pulled wa v e regime where the trav elling w av e sp eed is gov erned by the leading edge dynamics. The analysis that follows builds on earlier v ariational studies [ 21 , 25 ], but is reformulated here to yield an explicit condition. If D is contin uous near the origin, f ∈ C 1 near the origin, and D ( u ) f ( u )(1 − u ) β is b ounded, then the dominant contribution of the integral in the n umerator of F ( β ) as β ↗ 2 arises from a neigh b ourho od of u = 0 . In particular, if D (0) f ′ (0)  = 0 , then Z 1 0 D ( u ) f ( u ) u − β (1 − u ) β d u = D (0) f ′ (0) Z 1 0 u − β +1 d u + O (1) = D (0) f ′ (0) 2 − β + O (1) . Similarly , the Beta function satisfies B(2 − β , 2 + β ) = (2 − β ) − 1 − 11 / 6 + O (2 − β ) , and hence lim β ↗ 2 F ( β ) = 2 D (0) f ′ (0) , whic h reco vers the linear marginal stability b ound c ≥ c L := 2 p f ′ (0) D (0) . Since the minimal admissible wa v e sp eed is determined by c 2 = 2 sup β ∈ [0,2) F ( β ), w e see that nonlinear (pushed) wa v e sp eed selection o ccurs precisely when the supremum of F ( β ) is attained at an interior p oin t β < 2 . A sufficient lo cal condition for this is that F ( β ) increases as β decreases from β = 2 , namely d d β F ( β )     β =2 − < 0. (15) W e therefore study the sign of F ′ ( β ) as β ↗ 2 , assuming that F ′′ ( β ) ≤ 0 on [0, 2) so that F is conca ve. This criterion (Eq. (15) ) immediately explains several kno wn regimes. If f ( u ) ≥ 0 for all u ∈ [0, 1] and f ′ (0) D (0) = 0 , then F ( β ) ≥ 0 and the w av e speed is necessarily selected nonlinearly . This 7 includes P orous–Fisher-type mo dels with D ( u ) ∼ u m as u → 0 , as well as mo dels with Allee-t yp e reaction terms, where f ( u ) = u ( u − a )(1 − a ) is negative at low densities. When D (0) f ′ (0)  = 0 , the sign of F ′ (2 − ) determines whether nonlinearities in the diffusion or reaction terms shift the selected sp eed a wa y from c L . If F ′ ( β ) < 0 , then the sp eed is nonlinearly selected, as w e can find a low er b ound for the sp eed on the in terv al β ∈ (0, 2) , which exceeds c L . If we set R ( u ) through the relationship f ( u ) = uR ( u ) , w e can define N ( β ) = Z 1 0 D ( u ) f ( u ) u − β (1 − u ) β d u = Z 1 0 D ( u ) R ( u ) u 1 − β (1 − u ) β d u , with R (0)  = 0 . W e further assume that Z 1 0 D ( u ) R ( u ) − R (0) D (0) u d u < ∞ . In particular, this means that D ( u ) R ( u ) − D (0) R (0) = O ( u a ) for any a > 0 . In this case N ( β ) = D (0) R (0)B(2 − β , 1 + β ) + Z 1 0 D ( u ) R ( u ) − R (0) D (0) u β − 1 (1 − u ) β d u = D (0) R (0) 2 − β − 3 D (0) R (0) 2 + Z 1 0 D ( u ) R ( u ) − R (0) D (0) u β − 1 (1 − u ) β d u + O (2 − β ) . Using this expression for N ( β ) together with the asymptotic expansion of the Beta function, w e obtain F ( β ) = 2 D (0) R (0) | {z } F (2) − 2 3  1 2 D (0) R (0) − 3 Z 1 0 D ( u ) R ( u ) − R (0) D (0) u (1 − u ) 2 d u  | {z } F ′ (2) (2 − β ) + O ((2 − β ) 2 ) . Hence we obtain a sufficient condition for nonlinear sp eed selection ( i.e. F ′ (2) < 0 ) Z 1 0 D ( u ) R ( u ) − R (0) D (0) u (1 − u ) 2 d u > 1 6 D (0) R (0) , (16) where we note that R (0) = f ′ (0) . F or instance, if D ( u ) = u + δ and f ( u ) = u (1 − u ) , we hav e R ( u ) = 1 − u and D ( u ) R ( u ) − D (0) R (0) = u (1 − u − δ ) . Eq. (16) predicts a nonlinear contribution to the selected w av e sp eed when δ < 1 / 2 , agreeing with previous results [25]. The criterion given in Eq. (16) reflects a comp etition b et w een leading edge linearisation and nonlinear cont ributions in the wa ve in terior, and this mec hanism p ersists in multi-species systems under w eak or smo oth coupling—see Sec. 4. The remainder of this section fo cuses on applying the v ariational approach to a num b er of one-sp ecies models adapted from the literature to improv e the 8 estimated low er b ound on the trav elling wa v e sp eed. 3.2 A mo del of epidermal w ound healing. First, w e consider a generalised, non-dimensional equation which can b e used to study wound healing and other applications [ 6 , 28 , 29 ]. The gov erning reaction–diffusion equation takes the form ∂ t ρ = ∂ x ( ρ m ∂ x ρ ) + ρ (1 − ρ n ), (17) for p ositiv e integers m , n . This is an example of the well-studied Porous–Fisher equation, with a Ric hards’ gro wth term [ 12 , 30 ]. W e apply the v ariational principle to estimate the wa ve sp eed, using D ( u ) = u m , and f ( u ) = u (1 − u n ) . In this case, F ( β ) can b e calculated explicitly as F ( β ) = β B(2 − β , 2 + β ) " B( m − β + 2, β + 1) − B( m + n − β + 1, β + 1) # = 6 β (1 + β )Γ(2 − β ) " Γ( m − β + 2) Γ( m + 3) − Γ( m + n − β + 2) Γ( m + n + 3) # . (18) T aking the expression for F ( β ) from Eq. (18) and ev aluating the supremum in Eq. (11) giv es a new lo wer b ound on the tra velling wa ve sp eed, c , which cannot b e ev aluated explicitly , but is plotted n umerically in Fig. 1. 1 2 3 4 5 n 0 . 2 0 . 4 0 . 6 0 . 8 c m 1 2 3 4 5 1 2 3 4 5 m n 1 2 3 4 5 Figure 1: Plot of the numerically estimated minimum tra velling w av e sp eed (solid lines) and the minim um tra velling wa ve sp eed of Eq. (17) predicted using the v ariational principle with F ( β ) as given in Eq. (18) (dashed lines) as a function of n (left-hand side) and m (righ t-hand side). The numerically estimated trav elling w av e sp eed is obtained by tracing the p oin t X ( t ) suc h that u ( X ( t ), t ) = 0.1 . Sev eral k ey features are apparent in Fig. 1. First, across all parameter regimes shown, the v ariational principle pro vides a strict low er b ound, with the dashed curv es alw ays lying b elo w the n umerically observ ed minimum tra velling w av e sp eeds. The b ound is tightest when the diffusivity is closest to linear, i.e. , for small v alues of m , and when the reaction term is closest to logistic growth, i.e. , for small v alues of n . In these regimes, the wa v e profile b eha v es more similarly to the Fisher–KPP equation, for which the v ariational b ound is known to b e exact. F urthermore, as m increases, the 9 diffusion exhibits a stronger degeneracy , leading to a sharp er wa vefron t. In this limit, the v ariational b ound underestimates the tra velling wa ve sp eed more strongly as the impact of the diffusion term increases. A similar trend o ccurs when n increases: the reaction term b ecomes more strongly nonlinear, mo ving the system aw ay from the Fisher–KPP regime in whic h the v ariational b ound is sharp. In b oth cases, the gap b etw een the dashed and solid curv es reflects precisely the departure from the classical marginal-stability scenario that underpins the v ariational metho d. These trends motiv ate a closer examination of the sp ecial case n = 1 ( i.e. , the P orous–Fisher equation), for whic h the v ariational expression simplifies substantially and exact analytical results are a v ailable for several in teger v alues of m . When n = 1 , we use Eq. (11) to find F ( β ) = β B( m + 2 − β , 2 + β ) B(2 − β , 2 + β ) = β Γ(4) Γ(4 + m ) Γ(2 − β + m ) Γ(2 − β ) , (19) whic h can b e simplified for small integer v alues of m , as F ( β ) is a low-order p olynomial in β . F or instance, w e recov er the known v alues of the tra velling wa v e sp eed when m = 0 and m = 1 . F or m = 0 , Eq. (19) reduces to F ( β ) = β , and thus w e find c ≥ 2 by Eq. (11) (a.k.a. the Fisher–KPP tra velling w av e sp eed). F or m = 1 , using the fact that Γ( z + 1) = z Γ( z ) , w e see that c ≥ 1 / √ 2 , and for the case m = 2 , we find the explicit solution c ≥ q (10 + 7 √ 7) / 270 . W e highlight that this analysis also applies to the more general case where D ( u ) = u p , and f ( u ) = u q (1 − u ) , by defining m = p + q . The next example illustrates a contrasting scenario: a pushed-wa ve regime where the v ariational principle b ecomes exact. 3.3 A mo del with degenerate densit y-dep enden t diffusion and the Allee effect. In this example we consider a degenerate, parab olic reaction–diffusion equation, adapted from that studied in [ 17 ] to demonstrate how the v ariational approach can b e used to estimate the tra velling w av e sp eed when the Allee effect is present [31]. The go verning equation tak es the form ∂ t ρ = ∂ x (( αρ + ρ 2 ) ∂ x ρ ) + ρ (1 − ρ )( ρ − a ), (20) where D ( ρ ) = αρ + ρ 2 and f ( ρ ) = ρ (1 − ρ )( ρ − a ) for α > 0 and a ∈ [0, 0.5] . Employing the v ariational principle in Eq. (11), we find F ( β ) = β 120Γ(2 − β ) [6( α − a )Γ(4 − β ) + Γ(5 − β ) − 30 αa Γ(3 − β )] = β (2 − β ) 120 [(3 − β )(4 − β + 6( α − a )) − 30 α a ] , (21) where the recursion formula is used to simplify the Gamma functions. Emplo ying Eq. (21) in the expression for the minimum trav elling wa ve sp eed given in Eq. (11) , we numerically solve for the suprem um and find a very goo d prediction of the numerically observed tra velling w av e sp eed, which can b e observed in Fig. 2. It is notable that the accuracy of the v ariational b ound in this example is markedly higher than in the previous P orous–Fisher t yp e case, and indeed b ecomes exact for certain parameter regimes. This 10 0 1 2 α 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 c a 0 0.1 0.2 0.3 0.4 0.5 0 . 0 0 . 2 0 . 4 a α 0 0.5 1 1.5 2 2.5 Figure 2: Plot of the numerically estimated minimum tra velling w av e sp eed (solid lines) and the minim um tra velling wa ve sp eed of Eq. (20) deriv ed using the v ariational principle with F ( β ) found in Eq. (21) (dashed lines) as a function of α (left-hand side) and a (righ t-hand side). The numerically estimated trav elling w av e sp eed is obtained by tracing the p oin t X ( t ) such that u ( X ( t ), t ) = 0.1 . impro vemen t is not accidental: it reflects key structural differences introduced b y the Allee effect. When f ( u ) = ρ (1 − ρ )( ρ − a ) , the reaction term no longer exhibits logistic b eha viour and the classical marginal-stabilit y mechanism (where the wa ve sp eed is determined solely by the linearisation near u = 0 ) is replaced by a pushed-w av e regime in which the en tire in terior of the w av e profile contributes to sp eed selection. In contrast, for higher p o wers of the p orous–Fisher equation, no admissible test function can b e shown to attain equality in the v ariational principle. In biological terms, the Allee effect suppresses growth at lo w densities, shifting the leading-edge dynamics a wa y from the small u asymptotics that the v ariational principle relies on, but also, in turn, reduces the accuracy of the v ariational b ound. As a consequence, the v ariational metho d instead here captures the dominan t nonlinear contributions, pro ducing a near-exact estimate across the full parameter ranges of α and a , as shown in Fig. 2. 3.4 Fisher–Stefan mo del. In this example, the v ariational principle is adapted to study a w ell- kno wn mo ving b oundary problem, the Fisher–Stefan mo del, which is describ ed b y ∂ t ρ = ∂ 2 x ρ + ρ (1 − ρ ) , x < s ( t ) (22) ρ = 0, d s d t = − κ∂ x ρ , on x = s ( t ) , (23) where κ ≥ 0 is a parameter that controls the sp eed of the moving b oundary . By fixing the moving b oundary at z = 0 , w e can use the trav elling wa v e ansatz, z = x − s ( t ) = x − ct with ρ ( x , t ) = u ( z ) to obtain u ′′ + cu ′ + u (1 − u ) = 0 (24) u (0) = 0, u ′ (0) = − c κ . 11 As b efore, expressing v ( u ) = − d u/ d z as a function of u , and integrating Eq. (24) against a test function, φ ( u ) , we obtain Z 1 0  − 1 2 ϕ ( u ) v 2 ( u ) + cv ( u ) φ ( u )  d u + c 2 φ (0) 2 κ 2 = Z 1 0 u (1 − u ) φ ( u ) d u , where we used that v (0) = c/κ . Using the same argumen t as previously shown on the second order p olynomial in v , we obtain the b ound 1 2 c 2 ≥ R 1 0 u (1 − u ) φ ( u ) d u φ (0) /κ 2 + R 1 0 φ 2 ( u ) /ϕ ( u ) d u , with equality if v = cφ ( u ) ϕ ( u ) ⇐ ⇒ φ ( u ) ∝ exp  − Z u 0 c d q v ( q )  . Ev aluating at u = 0 , we also obtain the compatibilit y condition φ (0) = ϕ (0) /κ , whic h means that φ ( u ) = e − κu ˜ φ ( u ) with ˜ φ ′ (0) = 0 . Without loss of generality we assume φ (0) = ˜ φ (0) = 1 . The simplest test function we could employ that satisfies these b ounds and the compatibility condition is φ ( u ) = e − κu . Using this, w e find 1 2 c 2 ≥ R 1 0 u (1 − u ) e − κu d u 1 /κ 2 + R 1 0 e − κu /κ d u = κ − 2 + e − κ (2 + κ ) κ (2 + e − κ ) , (25) whic h giv es the correct b eha viour [32] as κ → 0 + : c ∼ κ √ 3 . Eq. (25) therefore provides a low er b ound for the trav elling wa ve sp eed that is v alid for all κ > 0 ; ho wev er, the quality of this b ound degenerates as κ → + ∞ (see Fig. 3, where it is notable that the b ound tends to one in the limit κ → ∞ while the sp eed should actually tend to tw o [33]). The degeneration of the b ound in the limit κ → ∞ is not just a shortcoming of the c hosen test function, but reflects a fundamental limitation of the v ariational form ulation in this regime. As κ increases, the influence of the Stefan b oundary condition weak ens and the Fisher–Stefan mo del approac hes the classical Fisher–KPP equation, for which the trav elling wa ve sp eed is selected b y linear marginal stability at the leading edge. The v ariational principle employ ed here enco des the mo ving b oundary condition through constraints on admissible test functions at u = 0 . In the limit κ → ∞ , these constraints degenerate, and the class of admissible test functions no longer captures the linearised dynamics gov erning sp eed selection. As a consequence, the resulting b ound saturates at a v alue strictly b elo w the true wa ve sp eed c = 2 . 12 − 10 − 5 0 z = x − ct 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 κ 0.25 2 5 10 20 50 0 5 10 15 20 κ 0 . 00 0 . 25 0 . 50 0 . 75 1 . 00 1 . 25 c Figure 3: (left) Shap e of trav elling w av e profiles for the Fisher-Stefan mo del giv en by Eqs. (22) and (23) dep ending on the parameter κ . (right) Plot of the numerically simulated minimum trav elling wa v e sp eed (squares) and analytically calculated tra velling w av e sp eed using the v ariational principle in Eq. (25) (blac k dashed lines) as a function of κ . The numerically estimated trav elling wa ve speed is obtained by tracing the p oin t X ( t ) such that u ( X ( t ), t ) = 0.1 . This observ ation highlights that, while the v ariational approach yields a rigorous and informative lo wer b ound for all κ > 0 , it is most effective when wa ve sp eed selection is gov erned by nonlinear or b oundary effects. Recov ering the Fisher–KPP sp eed in the large κ limit using these techniques would require a reformulation of the v ariational principle that explicitly incorp orates the linear leading edge b eha viour in this asymptotic regime, which lies b ey ond the scop e of the presen t work. 4 T w o-sp ecies models: a leading order v ariational appro ximation Although the full optimal control problem derived in Sec. 2 is exact, solving the asso ciated optimalit y system analytically is generally in tractable for mo dels inv olving tw o or more sp ecies. Instead, in this section, we demonstrate ho w we can exploit the optimal control approach to find a leading order v ariational approximation that captures the dominant wa ve sp eed selection mec hanism in a weak coupling regime b et ween the t wo species. The low er b ound for the minimum trav elling wa ve sp eed found using this v ariational approach impro ves the b ound found using linear stability analysis, and th us is particularly useful in sp ecific parameter regimes where further analytical progress is made p ossible. In particular, we restrict our attention to a particularly relev an t class of mo dels where ρ 2 denotes a densit y field which gets degraded by the first sp ecies, such that g ( ρ 1 , ρ 2 ) = − κρ 1 ρ 2 with ρ 2 (+ ∞ ) = ν . Solving for u 2 , we find u 2 ( u 1 ) = ν exp  − κ c Z u 1 0 q d q v ( q )  . (26) In this case, weak coupling corresp onds to small v alues of κ and ν , such that the secondary sp ecies remains uniformly small across the w av efron t. Indeed, Eq. (26) sho ws that u 2 ( u 1 ) = O ( ν ) uniformly on [0, 1] . W e therefore introduce the small parameter ϵ = κν /c , and seek an approximate solution 13 v alid as ϵ → 0 . W e start b y noting that for t wo-species systems, the optimal control v ∗ can b e found by solving the equations 0 = − ϕD 2 v 3 + cD φv 2 + κ c ω u 1 u 2 , (27) d u 2 d u 1 = − κ c u 1 u 2 v , (28) d ω d u 1 =  ϕD v 2 − cv φ  ∂ 2 D + ∂ 2 ( D f ) φ + κ c u 1 ω v , (29) with u 2 (0) = ν , u 2 (1) = 0 , and ∂ 2 = ∂ /∂ u 2 , D = D ( u 1 , u 2 ) and f = f ( u 1 , u 2 ) , for notational simplicit y . Motiv ated by the v ariational principle for single–sp ecies mo dels, we seek a solution of the form v = cφ/ϕD + δ v where δ v = O ( ϵ ) uniformly in [0, 1]. First, we observe that the adjoint enters the stationarit y condition (Eq. (27) ) only through the product u 2 ω . Differentiating this pro duct and using Eqs. (28) and (29) yields d d u 1 ( u 2 ω ) = u 2  ϕD v 2 − cv φ  ∂ 2 D + u 2 ∂ 2 ( D f ) φ , and hence u 2 ( u 1 ) ω ( u 1 ) = − Z 1 u 1  u 2  φD v 2 − cv φ  ∂ 2 D + u 2 ∂ 2 ( D f ) φ  d q . Since the secondary species satisfies 0 ≤ u 2 ( u 1 ) ≤ ν uniformly on [0, 1] , and admissible test functions ob ey φ ( u 1 ) ≲ u − β 1 with β ∈ (0, 2) , the integral representation for the adjoint con tribution u 2 ω remains uniformly b ounded on [0, 1] . As a consequence, u 2 ω = O ( ν ) uniformly , and the adjoin t term in the stationarit y condition (Eq. (27)) enters only at order O ( κν /c ) = O ( ϵ ) . Under mild regularity assumptions on the reaction and diffusion functions, and standard ad- missibilit y conditions on the test function, we show in App endix A that the ansatz v = cφ/ϕ + δ v satisfies the full optimality condition (Eq. (27) ) up to order O ( ϵ ) . As such, at leading order in the w eak coupling parameter ϵ , the optimal control form ulation coincides with that of the one-sp ecies v ariational principle, given by 1 2 c 2 ≥ R 1 0 D ( u 1 , u 2 ( u 1 )) f ( u 1 , u 2 ( u 1 )) φ ( u 1 ) d u 1 R 1 0 φ ( u 1 ) 2 /ϕ ( u 1 ) d u 1 , with ϕ = − φ ′ , (30) where u 2 ( u 1 ) is given b y Eq. (26) and equalit y is only achiev ed when v ∗ = cφ ( u 1 ) D ( u 1 , u 2 ) ϕ ( u 1 ) . (31) This leading order reduction provides a practical v ariational b ound for m ultisp ecies systems. In the following subsections, we apply this framework to biologically motiv ated mo dels and demonstrate that it captures nonlinear wa v e sp eed selection b ey ond linear theory . 14 4.1 Mo dels of cell in v asion into extracellular matrix. W e now consider a class of mo dels describing collectiv e cell migration into the extracellular matrix (ECM) which take the following form ∂ t ρ 1 = ∂ x ((1 − ρ 2 ) ∂ x ρ 1 ) + f ( ρ 1 , ρ 2 ) , (32) ∂ t ρ 2 = − κρ 1 ρ 2 , (33) where f (0, ρ 2 ) = 0 and f (1, ρ 2 ) = 0. The b oundary conditions are given by    ρ 1 ( x , t ) → 1, ρ 2 ( x , t ) → 0, as x → −∞ , ρ 1 ( x , t ) → 0, ρ 2 ( x , t ) → ν , as x → + ∞ , for 0 ≤ ν ≤ 1. W e consider tw o different reaction terms. Firstly , f C ( ρ 1 , ρ 2 ) = f C ( ρ 1 ) = ρ 1 (1 − ρ 1 ) as describ ed in [ 11 ]. Secondly , w e consider f B ( ρ 1 , ρ 2 ) = ρ 1 (1 − ρ 1 − ρ 2 ) as detailed in [ 5 , 13 ]. F or b oth mo dels, n umerical simulations show that the system exhibits trav elling w av es that attain the linearly predicted w av e sp eed ( 2 √ 1 − ν when f = f C , and 2(1 − ν ) when f = f B ) for small v alues of the ECM degradation rate, κ ≤ κ ∗ ( ν ) , only , where κ ∗ ( ν ) is the critical ECM degradation rate such that the linearly predicted and numerically observed trav elling wa v e sp eeds coincide for κ ≤ κ ∗ ( ν ) . F or sufficien tly large ECM degradation rates κ > κ ∗ ( ν ) , this mo del exhibits nonlinear wa ve sp eed selection. The associated numerically estimated wa ve sp eed sho w significant deviations from the linearly selected tra velling wa ve sp eed, which underestimates the in v asion sp eed. F urthermore, the critical v alue of the degradation rate, κ ∗ ( ν ) , seems to follow a monotonically increasing relationship with ν as ν → 1 [ 9 ]. W e now wan t to apply the leading order v ariational approximation from Sec. 4. Employing the same family of test functions as in Sec. 3, ( φ β ( u 1 ) = ((1 − u 1 ) /u 1 ) β for β ∈ [0, 2) ), we ha ve 1 2 c 2 ≥ sup β ∈ [0,2) G ( β ), G ( β ) = β M ( β ) B(2 − β , 2 + β ) , (34) where M ( β ) = Z 1 0 D ( u 1 , u 2 ( u 1 )) f ( u 1 , u 2 ( u 1 ))  1 − u 1 u 1  β d u 1 . (35) F or D ( u 1 , u 2 ( u 1 )) = (1 − u 2 ( u 1 )) and f C ( u 1 ) = u 1 (1 − u 1 ) , we find M C ( β ) = Z 1 0 u 1 − β (1 − u ) 1+ β ν 1 − ν (1 − u ) κβ /c 2 +1 d u . (36) Ho wev er, for f B ( u 1 , u 2 ) = u 1 (1 − u 1 − u 2 ) with γ = ν / (1 − ν ) we hav e M B ( β ) = Z 1 0 (1 − u 2 ( u 1 ))(1 − u 1 − u 2 ( u 1 )) u 1 − β 1 (1 − u 1 ) β d u 1 = Z 1 0 u 1 − β 1 (1 − u 1 ) β − u 2 − β 1 (1 − u 1 ) β − γ (1 − u 1 ) 1 − β (1 − u 1 ) β + η (1 + γ (1 − u 1 ) η ) 2 . (37) 15 10 ° 3 10 ° 1 10 1 10 3 10 5 ∑ 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 c ∫ 0.0 0.2 0.4 0.6 0.8 1.0 10 ° 3 10 ° 1 10 1 10 3 10 5 ∑ 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 c ∫ 0.0 0.2 0.4 0.6 0.8 1.0 A B Figure 4: Plot of n umerically simulated minimum trav elling wa ve sp eed (solid lines) and analytically calculated tra velling wa ve sp eed of Eqs. (32) and (33) using linear theory (dotted lines) and the v ariational principle (dashed lines) with M C ( β ) as calculated in Eq. (36) (left-hand side) and M B ( β ) as calculated in Eq. (37) (righ t-hand side) for v arious initial ECM densities ν ≥ 0 . The n umerically estimated trav elling wa ve sp eed is obtained by tracing the p oint X ( t ) such that u 1 ( X ( t ), t ) = 0.1 . Substituting Eqs. (36) and (37) in to Eq. (30) pro vides a low er b ound on the tra velling wa ve sp eed that, in each case, can b e solved for n umerically . Fig. 4A shows the result of n umerically solving Eq. (30) sub ject to Eq. (36) (a.k.a. when f ( u 1 , u 2 ) = f C ( u 1 ) ). Similarly , Fig. 4B sho ws the result of n umerically solving Eq. (30) sub ject to Eq. (37) (a.k.a. when f ( u 1 , u 2 ) = f B ( u 1 , u 2 ) ) for a low er b ound on the trav elling wa v e sp eed, c , whic h shows a v ast improv ement when compared to the linearly predicted tra velling wa v e sp eeds. T o the b est of our knowledge, the linear tra velling wa v e sp eed is the curren t b est estimate of the sp eed of inv asion in the literature for b oth of the mo dels presen ted in this section, so the prediction deriv ed in this w ork using the v ariational approach drastically impro ves the low er b ound on the wa ve sp eed in b oth cases. 4.2 A mo del of cell in v asion with cell differen tiation. W e no w consider a t wo–species model of cell inv asion first introduced in [ 3 ], which has recen tly attracted renewed interest due to its transition from linear to nonlinear wa ve sp eed selection [ 19 ]. The mo del describ es a p opulation of inv asive cells, c( x , t ) , that migrate via diffusion and proliferate logistically , while differentiating at rate λ > 0 into a p opulation of non–motile, non–proliferative cells, n ( x , t ) —see Fig. 5 for n umerical solutions. The go verning equations are ∂ t c = ∂ 2 x c + c(1 − c) − λ c(1 − K n ) , (38) ∂ t n = λ c(1 − K n ) , (39) for K > 0 mo dulating the differen tiation rate. The steady states of this system are (c, n ) = (0, n ∗ ) and (1, 1 /K ) , for an y 0 ≤ n ∗ ≤ 1 /K , where w e choose n ∗ = 0 as this corresp onds to the more natural case of having no differentiated cells initially [ 3 ]. This mo del can b e reduced to a form similar to that 16 in Eqs. (32) and (33) by setting ρ 2 = 1 − K n , ρ 1 = c . This gives ∂ t ρ 1 = ∂ 2 x ρ 1 + ρ 1 (1 − ρ 1 − λρ 2 ) , (40) ∂ t ρ 2 = − κρ 1 ρ 2 , (41) where κ = λK . W e assume that λK ≪ 1 so that we can use the leading order v ariational principle w e deriv ed in Sec. 4. The b oundary conditions of the system are given by    ρ 1 ( x , t ) → 1 ρ 2 ( x , t ) → 0, as x → −∞ , ρ 1 ( x , t ) → 0 ρ 2 ( x , t ) → 1, as x → + ∞ . Setting ρ 1 ( x , t ) = u 1 ( x − ct ) and ρ 2 ( x , t ) = u 2 ( x − ct ) , the linear marginal stabilit y condition in this case giv es c ≥ c L = 2 √ 1 − λ whic h only applies for λ ≤ 1 . Emplo ying the leading order v ariational principle using Eq. (35), we calculate M ( β ) = Z 1 0 u 1 − β 1 (1 − u 1 ) β h (1 − u 1 ) − λ (1 − u 1 ) β κ/c 2 i d u 1 = B(2 − β , 2 + β ) − λB  2 − β , 1 + β + β κ c 2  . Hence, Eq. (34) giv es G ( β ) = β   1 − 6 λ Γ  1 + β + β κ c 2  Γ  3 + β κ c 2  Γ(2 + β )   . (42) W e note that G (2) = 2(1 − λ ) , giving c 2 ≥ 4(1 − λ ) . This is in agreement with the trav elling w av e sp eed predicted by linear theory . In particular, our approach recov ers, and in fact improv es up on, the marginal b ound obtained from linear stability analysis. Next, w e use Eq. (42) to obtain the sharpest b ound on the w a ve sp eed in the form c 2 / 2 ≥ sup β ∈ (0,2) G ( β ) . In Fig. 5 we also indicate the parameter regimes in which the weak coupling appro xi- mation breaks down, namely when λκ/c ∼ 1 . Overall, we observe that the approximate v ariational principle provides an accurate prediction of the wa ve sp eed. Agreemen t is excellent within the weak coupling regime, and although deviations b ecome visible as this regime is exited, the resulting error remains comparatively small. It is particularly in teresting that the approximation remains accurate even for large v alues of K . This motiv ates an in vestigation of the cell differen tiation mo del in the regime where κ = λK is large but λ remains small. In this case the adjoin t v ariable ω satisfies ω ( u 1 ) u 2 ( u 1 ) = λ Z 1 u 1 q φ ( q ) u 2 ( q ) d q . It is easy to sho w that when κ → + ∞ and λ ≪ 1 the ansatz v ∗ = cφ/ϕ is the leading order 17 0 100 z = x − ct 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 λ = 0 . 5 , K = 3 c n 0 100 z = x − ct λ = 1 . 5 , K = 3 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 λ 0 . 0 0 . 5 1 . 0 1 . 5 2 . 0 c c = 2 √ 1 − λ K = 0 . 5 K = 1 K = 3 K = 10 K = 50 Figure 5: (left) Numerical simulations of the time-dep enden t cell differentiation model (Eqs. (38) and (39) ) exhibiting trav elling w av e solutions. Solutions plotted at t = 50, 75, 100 for different parameters. (righ t) W av e sp eed c as a function of λ for v arious v alues of K . Solid lines corresp ond to predictions obtained from the appro ximate v ariational principle, while square markers indicate numerical solutions of the full PDE system. Dashed vertical lines mark the threshold λκ/c = 1 , b ey ond which the weak coupling approximation is no longer v alid. The blac k curv e shows the linear theory prediction c = 2 √ 1 − λ . solution to the optimal control formulation. Indeed, using this c hoice of v to solv e for u 2 , with φ ( u 1 ) = ((1 − u 1 ) /u 1 ) β , and expanding for large κ , we obtain ω ( u 1 ) u 2 ( u 1 ) = λ Z 1 u 1 q 1 − β (1 − q ) β + κβ /c 2 d q ∼ u 1 − β 1 (1 − u 1 ) κβ /c 2 β + κβ /c 2 . In particular, as κ → + ∞ , we obtain κ c ω u 1 u 2 = O ( cλ ) , (43) so the weak coupling assumption is in fact satisfied if λc is small, which explains the go od agreement observ ed in Fig. 5. 4.3 Extension to mo dels with taxis terms. Benguria et al. ha ve, more recently , considered the addition of advection terms in the gov erning equations and their impact on the resulting trav elling w av e sp eed [ 23 ]. As is natural, we therefore presen t a simple extension of the work presented here that is applicable to p opulations undergoing taxis, or following a gradient of external signals, for example, chemicals or substrates— e.g. , chemotaxis, haptotaxis, electrotaxis. The mo dels of interest will take the form ∂ t ρ 1 = ∂ x ( D ( ρ 1 , ρ 2 ) ∂ x ρ 1 − χ ( ρ 1 , ρ 2 ) ∂ x ρ 2 ) + f ( ρ 1 , ρ 2 ) , (44) ∂ t ρ 2 = g ( ρ 1 , ρ 2 ) , (45) where the b oundary conditions of the system are given b y    ρ 1 ( x , t ) → 0 ρ 2 ( x , t ) → µ , as x → −∞ , ρ 1 ( x , t ) → 1 ρ 2 ( x , t ) → ν , as x → + ∞ . 18 If w e ha ve that D ( ρ 1 , ρ 2 ) ≥ 0 is the density-dependent diffusivity , χ ( ρ 1 , ρ 2 ) ≥ 0 is the tactic sensitivit y , and f ( ρ 1 , ρ 2 ), g ( ρ 1 , ρ 2 ) are reaction terms, w e can then lo ok for a monotonically decreasing function of z = x − ct , defining ρ 1 ( x , t ) = u 1 ( x − ct ) and ρ 2 ( x , t ) = u 2 ( x − ct ) , such that    u 1 ( z ) → 0 u 2 ( z ) → µ as z → −∞ , u 1 ( z ) → 1 u 2 ( z ) → ν as z → ∞ . If w e consider a w eak coupling with the external cue, we can again follo w the v ariational principal approac h detailed earlier in Sec. 4. Namely , setting d u 1 / d z = − v , w e ha ve d u 1 d z = − v , d d z  D ( u 1 , u 2 ) v + χ ( u 1 , u 2 ) d u 2 d z  = f ( u 1 , u 2 ) − cv , d u 2 d z = − 1 c g ( u 1 , u 2 ), Replacing all instances of D ( u 1 , u 2 ) v b y D ( u 1 , u 2 ) v + χ ( u 1 , u 2 ) d u 2 / d z in the calculations, we can follo w our previous w orking to sho w an approximate v ariational principle for c hemotactic style models. Ultimately , we find 1 2 c 2 ≥ R 1 0 f ( u 1 , u 2 ( u 1 )) φ ( u 1 )  D ( u 1 , u 2 ) − χ ( u 1 , u 2 ) d u 2 d u 1  d u 1 R 1 0 φ ( u 1 ) 2 /ϕ ( u 1 ) d u 1 , with ϕ = − φ ′ , as long as d u 2 d u 1 = ϕ ( u 1 )  D ( u 1 , u 2 ) − χ ( u 1 , u 2 ) d u 2 d u 1  g ( u 1 , u 2 ) c 2 φ ( u 1 ) . (46) Due to lac k of examples in the literature, this metho d is not applied to any mo dels in this work but instead detailed to demonstrate the utility and wider applicability of the approach to these mo dels, where they exist. 5 Discussion In this article, w e ha ve developed an optimal control framework for c haracterising the minimal tra velling w av e sp eed in systems of reaction–diffusion equations with nonlinear diffusion or reaction terms admitting monotonic fron ts. This approach yields a v ariational characterisation of the minimal sp eed that reco vers classical results for single–sp ecies equations [ 21 , 23 – 25 ] and, in a weak coupling regime, extends them naturally to coupled systems. Within this framework, w e also deriv e a criterion for the onset of nonlinear w av e sp eed selection, and apply the asso ciated v ariational principle to sev eral biologically motiv ated mo dels to elucidate parameter dep endence. F rom an optimal control p erspective, the classical single–sp ecies v ariational principle emerges as a sp ecial, explicitly solv able case of a more general control problem: the tra velling wa ve profile 19 serv es as the state v ariable, the lo cal flux as the control, and the w av e sp eed arises as the optimal v alue satisfying a stationarity condition. This in terpretation clarifies the structure underlying existing v ariational argumen ts and provides a systematic route to generalising wa v e sp eed selection theory to m ulti-sp ecies systems, where closed-form v ariational formulations are t ypically una v ailable. Despite the b enefits of the v ariational characterisation, further research is needed to address the limitations of the metho ds. The analysis assumes the existence of monotonic trav elling wa v es with constan t sp eed, enabling the reduction of the PDEs to a system of ordinary differential equations in the tra velling wa ve co ordinate. Moreo ver, the multi-species v ariational approximation relies on sufficien tly w eak coupling so that the adjoint contributions in the optimalit y system may b e treated as higher-order corrections only . In regimes with strong coupling, sharp internal la yers, or strong dep endence of the diffusivit y of the motile sp ecies on the additional sp ecies, the resulting low er b ounds are not exp ected to b e sharp and may only capture qualitative parameter dep endence. Nonetheless, the framew ork presented here provides improv ed estimates on the lo wer b ound for the trav elling wa v e sp eed together with a unified and extensible foundation for studying nonlinear wa ve sp eed selection in general in systems of reaction–diffusion equations. A Uniform b oundedness of the adjoin t contribution In order to show that if v = cφ/D ϕ + δ v , where δ v = O ( ϵ ) , then v satisfies Eq. (27) uniformly in u 1 ∈ [0, 1] up to order O ( ϵ ) , we make the following mild assumptions on the reaction and diffusivit y of the first sp ecies, and on the test function. W e assume that there exist constants a 1 , a 2 , a 3 , a 4 > 0 suc h that max 0 ≤ u 1 ≤ 1, 0 ≤ u 2 ≤ ν D , | ∂ 2 D | , | ∂ 2 ( D f ) | ≤ a 1 , D ≥ a 2 > 0; (47) and a 3 (1 − u 1 ) β ≤ φ ≤ a 4 u − β , a 4 (1 − u 1 ) β − 1 ≤ ϕ ≤ a 4 u − β − 1 , for β ∈ (0, 2) . Consider the final term, inv olving ω , in Eq. (27) . Using the b oundedness assumptions, we find that | u 1 u 2 ω | ≤ C  ν c 2 u 1 Z 1 u 1 q 2 (1 − q ) 2 d q + ν u 1 Z 1 u 1 u − β  + O ( ϵ ) . And from here w e obtain sup β ∈ (0,2) max u 1 ∈ [0,1] | u 1 u 2 ω | ≤ ν C + O ( ϵ ) . where C = C ( a 1 , a 2 , a 3 , a 4 , c ) > 0 . Thus v ∗ = cφ/D ϕ solv es the optimal control problem, uniformly in u 1 ∈ [0, 1] , up to order O ( ϵ ) . 20 A c knowledgmen ts The authors would like to thank Mat Simpson and Philip Maini for helpful discussions. F or the purp ose of Op en A ccess, the authors hav e applied a CC BY public copyrigh t licence to any Author A ccepted Manuscript (AAM) version arising from this submission. R.M.C. would like to thank the Engineering and Ph ysical Sciences Research Council (EP/T517811/1) for funding. C.F. ackno wledges supp ort from a Ho ok e Researc h F ellowship. R.E.B. ackno wledges supp ort of a gran t from the Simons F oundation (MP-SIP-00001828). References 1. J. D. Murray . Mathematic al Biolo gy . Springer, German y , 1993. 2. M. J. Simpson, D. C. Zhang, M. Mariani, K. A. Landman, and D. F. Newgreen. Cell proliferation driv es neural crest cell in v asion of the intestine. Developmental Biolo gy , 302(2):553–568, 2007. 3. A. J. T rew enack and K. A. Landman. A trav eling wa v e mo del for inv asion by precursor and differen tiated cells. Bul letin of Mathematic al Biolo gy , 71:291–317, 2009. 4. P . Gerlee and S. Nelander. T ra velling w av e analysis of a mathematical mo del of glioblastoma gro wth. Mathematic al Bioscienc es , 276:75–81, 2016. 5. A. P . Browning, P . Haridas, and M. J. Simpson. A Ba yesian sequential learning framew ork to parameterise contin uum mo dels of melanoma inv asion into h uman skin. Bul letin of Mathematic al Biolo gy , 81(3):676–698, 2019. 6. C. F alcó, D. J. Cohen, J. A. Carrillo, and R. E. Bak er. Quan tifying tissue gro wth, shap e and collision via contin uum mo dels and Bay esian inference. Journal of the R oyal So ciety Interfac e , 20(204):20230184, 2023. 7. R. A. Fisher. The wa ve of adv ance of adv antageous genes. Annals of Eugenics , 7(4):355–369, 1937. 8. A. N. Kolmogoro v, I. Petro vskii, and N. S. Piskunov. A study of the equation of diffusion with increase in the quantit y of matter, and its application to a biological problem. Mosc ow University Biolo gic al Scienc es Bul letin , 1(6):1–25, 1937. 9. R. M. Crossley , P . K. Maini, T. Lorenzi, and R. E. Bak er. T rav eling wa ves in a coarse-grained mo del of volume-filling cell inv asion: Sim ulations and comparisons. Studies in Applie d Mathematics , 151(4):1471–1497, 2023. 10. E. Bouin and V. Calvez. T rav elling wa v es for the cane toads equation with bounded traits. Nonline arity , 27(9):2233, 2014. 11. C. Colson, F. Sánchez-Garduño, H. M. Byrne, P . K. Maini, and T. Lorenzi. T rav elling-w av e analysis of a mo del of tumour in v asion with degenerate, cross-dep enden t diffusion. Pr o c e e dings of the R oyal So ciety A , 477(2256):20210593, 2021. 21 12. A. De Pablo and A. Sánc hez. T rav elling wa ve b eha viour for a p orous–Fisher equation. Eur op e an Journal of Applie d Mathematics , 9(3):285–304, 1998. 13. M. El-Hachem, S. W. McCue, and M. J. Simpson. T rav elling wa ve analysis of cellular inv asion in to surrounding tissues. Physic a D: Nonline ar Phenomena , 428:133026, 2021. 14. C. F alcó, R. M. Crossley , and R. E. Bak er. T rav elling wa ves in a minimal go-or-grow mo del of cell inv asion. Applie d Mathematics L etters , 158:109209, 2024. 15. K.-P . Hadeler and F. Rothe. T rav elling fronts in nonlinear diffusion equations. Journal of Mathematic al Biolo gy , 2(3):251–263, 1975. 16. B. P . Marc han t, J. Norbury , and J. A. Sherratt. T ra velling w av e solutions to a haptotaxis- dominated mo del of malignant in v asion. Nonline arity , 14(6):1653, 2001. 17. F. Sánchez-Garduño and P . K. Maini. T rav elling wa v e phenomena in some degenerate reaction– diffusion equations. Journal of Differ ential Equations , 117(2):281–319, 1995. 18. T. L. Stepien, E. M. Rutter, and Y. Kuang. T ra veling w av es of a go-or-grow mo del of glioma gro wth. SIAM Journal on Applie d Mathematics , 78(3):1778–1801, 2018. 19. Y. Y ue and C. Ou. T rav eling w av efronts to a mo del of precursor and differentiated cells: Impacting parameter-structure transition from monostable to bistable, and from monotone to non-monotone. Journal of Dynamics and Differ ential Equations , 37(4):3043–3078, 2025. 20. W. V an Saarlo os. F ron t propagation into unstable states. I I. Linear versus nonlinear marginal stabilit y and rate of con vergence. Physic al R eview A , 39(12):6367, 1989. 21. R. D. Benguria and M. C. Depassier. Sp eed of fronts of the reaction–diffusion equation. Physic al R eview L etters , 77(6):1171, 1996. 22. R. D. Benguria and M. C. Depassier. V ariational characterization of the sp eed of propagation of fron ts for the nonlinear diffusion equation. Communic ations in Mathematic al Physics , 175:221–227, 1996. 23. R. D. Benguria, M. C. Depassier, and S. Rica. V ariational estimates for the sp eed propagation of fron ts in a nonlinear diffusiv e Fisher equation. Chaos, Solitons & F r actals , 164:112668, 2022. 24. R. D. Benguria and M. C. Depassier. V alidit y of the linear sp eed selection mec hanism for fronts of the nonlinear diffusion equation. Physic al R eview L etters , 73(16):2272, 1994. 25. B. M. Stok es, T. Rogers, and R. James. Sp eed and shap e of p opulation fronts with densit y- dep enden t diffusion. Bul letin of Mathematic al Biolo gy , 86(12):1–18, 2024. 26. S. Lenhart and J. T. W orkman. Optimal Contr ol Applie d to Biolo gic al Mo dels . Chapman and Hall/CR C, 2007. 22 27. R. Benguria, M. Depassier, and V. Méndez. Minimal sp eed of fronts of reaction–conv ection– diffusion equations. Physic al R eview E , 69(3):031106, 2004. 28. Y. Liu, K. Suh, P . K. Maini, D. J. Cohen, and R. E. Baker. P arameter iden tifiability and mo del selection for partial differential equation mo dels of cell inv asion. Journal of the R oyal So ciety Interfac e , 21(212):20230607, 2024. 29. J. A. Sherratt and J. D. Murray . Models of epidermal wound healing. Pr o c e e dings of the R oyal So ciety of L ondon. Series B: Biolo gic al Scienc es , 241(1300):29–36, 1990. 30. T. P . Witelski. Merging trav eling wa v es for the p orous–Fisher’s equation. Applie d Mathematics L etters , 8(4):57–62, 1995. 31. P . A. Stephens, W. J. Sutherland, and R. P . F reckleton. What is the Allee effect? Oikos , pages 185–190, 1999. 32. M. El-Hachem, S. W. McCue, and M. J. Simpson. In v ading and receding sharp-fronted tra velling w av es. Bul letin of Mathematic al Biolo gy , 83(4):35, 2021. 33. S. W. McCue, M. El-Hac hem, and M. J. Simpson. Exact sharp-fron ted tra velling w av e solutions of the Fisher–KPP equation. Applie d Mathematics L etters , 114:106918, 2021. 23