Distinguishing Between Long-Transient and Asymptotic States in a Biological Aggregation Model

Aggregations are emergent features common to many biological systems. Mathematical models to understand their emergence are consequently widespread, with the aggregation–diffusion equation being a prime example. Here we study the aggregation–diffusion equation with linear diffusion in one spatial dimension. This equation is known to support solutions that involve both single and multiple aggregations. However, numerical evidence suggests that the latter, which we term ‘multi-peaked solutions’ may often be long-transient solutions rather than asymptotic steady states. We develop a novel technique for distinguishing between long transients and asymptotic steady states via an energy minimisation approach. The technique involves first approximating our study equation using a limiting process and a moment closure procedure. We then analyse local minimum energy states of this approximate system, hypothesising that these will correspond to asymptotic patterns in the aggregation–diffusion equation. Finally, we verify our hypotheses through numerical investigation, showing that our approximate analytic technique gives good predictions as to whether a state is asymptotic or transient. Overall, we find that almost all twin-peaked, and by extension multi-peaked, solutions are transient, except for some very special cases. We demonstrate numerically that these transients can be arbitrarily long-lived, depending on the parameters of the system.


Introduction
Aggregation phenomena are widespread in biology, from cell aggregations (Budrene and Berg, 1995) to the swarming (Roussi, 2020), schooling (Makris et al, 2009), flocking (Clark and Mangel, 1984), and herding (Bond et al, 2019) of animals.When modelled from a continuum perspective (as opposed to via interacting particles), the principal tools take the form of partial differential equations with non-local advection, sometimes combined with a diffusive term (Topaz et al, 2006).Indeed, such equations are often called aggregation equations (Laurent, 2007), highlighting their importance in modelling aggregations, or aggregation-diffusion equations (Carrillo et al, 2019) if there is a diffusion term.
As well as modelling aggregated groups of organisms, such equations have also been used to model aggregation-like phenomena elsewhere, such as animal home ranges and territories (Briscoe et al, 2002;Potts and Lewis, 2016) and consensus convergence in opinion dynamics (Garnier et al, 2017).This very broad range of applications, together with the mathematical complexity in dealing with nonlinear nonlocal partial differential equations (PDEs), has led to a great amount of interest from applied mathematicians in understanding the properties of these PDEs (Painter et al, 2023).
Of particular interest from a biological perspective are the pattern formation properties of aggregation-diffusion equations, since these can reveal the necessary processes required for observed patterns to emerge.Many traditional techniques for analysing pattern formation, such as linear stability analysis and weakly nonlinear analysis, focus on the onset of patterns from small perturbations of a non-patterned (i.e.spatially homogeneous) state.However, patterns observed in actual biological systems will often be far from the non-patterned state, and not necessarily emerge from small perturbations of spatially homogeneous configurations (Krause et al, 2020;Veerman et al, 2021).
Sometimes observed patterns will be asymptotic steady states or other types of attractors.But frequently biological systems will be observed in transient states (Hastings et al, 2018;Morozov et al, 2020).These transient states may persist for a very long time, sometimes so long that they are hard to distinguish from asymptotic states.Moreover, as well as transients being difficult to decipher from observations of biological systems, they can also be tricky to determine from numerical solutions of a PDE model.Therefore analytic techniques are required to guide those engaging in numerical analysis of PDEs as to whether the solution they are observing is likely to be a long transient or an asymptotic state.
Our aim here is to provide such analytic techniques for a class of 1D aggregationdiffusion equations of the following form where K is a non-negative averaging kernel, symmetric about 0, with ∥K∥ L ∞ = 1, and is a convolution, where Ω is the spatial domain of definition.Here, D and γ are constants, and Ω is the circle given by interval [−L, L] with periodic boundary conditions imposed.
Our approach is not exact, in the sense that we approximate our study PDE first through the limit as D/γ → 0, then via a moment closure assumption.However, this approximation allows us to analyse the associated energy functional, finding explicit mathematical expressions for local energy minima.Our conjecture is that local energy minima of the approximate system are qualitatively similar to the asymptotic patterns observed the aggregation-diffusion equation we are studying, but any states that do not represent local energy minima of the approximate system are transient states.We then test this numerically in some specific cases.
Of particular interest is the question of whether multi-peaked solutions are asymptotic steady states or long transients, which is the question that originally motivated this work.Various numerical studies of Equation ( 1), and similar equations, report multi-peaked solutions (Armstrong et al, 2006;Buttenschön and Hillen, 2020;Carrillo et al, 2019;Daneri et al, 2022).However, merging and decaying of peaks have also been observed.Furthermore, analytic investigations into chemotaxis equations, which have some similarities with aggregation equations, have demonstrated that multi-peaked solutions can often be long transients (Potapov and Hillen, 2005).
This work demonstrates that, except for the very specific case where peaks are of identical heights and evenly-spaced, any two-peaked solutions will eventually evolve into a solution with at most one peak, as the smaller peak decays to zero.The time it takes for the smaller peak to decay grows rapidly with the start height of the smaller peak, eventually tending to infinity as the difference in start heights between the two peaks tends to zero.We show that a key parameter governing the speed of this decay is the diffusion constant D, with higher diffusion constants leading to faster decays.We conjecture that, as D → 0, the time to decay tends to infinity, meaning that two-peaked solutions become stable.
Finally, we investigate the effect of incorporating logistic growth of the population into our model.The motivation for this is that, in situations where transient solutions exists for a long time, it is no longer biologically reasonable to assume that we are working in situations where births and deaths are negligible.We show that, for a given set of parameters and initial condition, there is a critical net reproduction rate, below which the smaller peak will decay and above which it will persist.

Methodological approach
Our study is motivated by an observation.Often, when simulating Equation (1), multiple aggregations may form and persist for a very long time.This can give the appearance of multi-peaked asymptotic stable states.For example, Figure 1 shows a numerical solution where two peaks have formed by time t = 1.These appear stable on timescales up to two orders of magnitude longer than the time they took to form: even by time t = 100, the solution has not changed very much (Figure 1a).However, if we keep running the simulation, we see one of the peaks decay and the other slowly  19) with δ = 0.1.
swallow up the former's mass.The question then arises whether multi-peaked solutions to Equation (1) are ever actually stable, or are they always just long transients?
To answer this question, our approach will not be to analyse Equation (1) directly, but rather take two approximations, which enable us to perform analytic calculations.First, we assume that γ ≫ D. Second, we make the following moment closure assumption where is the second moment of K.This leads to the following approximate version of Equation ( 1) Note that Equations ( 1) and ( 5) both preserve mass when solved with periodic boundary conditions (i.e.u( for all t > 0. Our tactic will be to search for minimum energy solutions to Equation (5) using the following energy functional In particular, we are interested in examining critical points of E[u], so calculate Here, the second and fourth equalities use integration by parts, together with the periodic boundary conditions.If we assume that there exist non-negative solutions to Equation ( 5) then the final expression in Equation ( 9) is non-positive, so that E[u] is non-increasing.Whilst we do not currently have a proof of the non-negativity of u, we note that that all our numerics suggest that non-negativity is preserved over time, that non-negativity results exist for Equation (1) for a variety of different kernels K (Carrillo et al, 2019;Giunta et al, 2022a;Jüngel et al, 2022), and so conjecture these might be transferable to the situation of Equation ( 5) with some effort.Equation ( 9) shows that critical points, u * (x), of the energy functional occur when which means that, on any connected subset of [−L, L], either u * (x) = 0 or for constants A, B, and C. Numerics suggest that Equation (1) tends towards a solution containing one or many aggregations, interspersed by constant sections close or near to zero (e.g. Figure 1).We want to construct differentiable solutions that have this type of qualitative appearance, yet also correspond to critical points of E[u].These can be constructed piecewise from Equation (11).For example, as long as πσ < √ 2L, a single-peaked solution can be given as follows where ϵ ∈ 0, p 2L and c ϵ are constants.One can also construct multi-peaked solutions in a similar way (which we will do later in the case of two peaks).Notice that such solutions are continuously differentiable, i.e. u * ∈ C 1 ([−L, L]), but not necessarily twice differentiable, so need to be understood in a weak sense (Evans, 2022).
By Equation ( 7), a direct calculation gives so that the only free parameter in Equation ( 12) is ϵ.Since the energy, E[u], is nonincreasing over time, the question arises as to which value of ϵ minimises E[u] across the set of all functions of the form in Equation ( 12).Our approach is to derive such minima, both in the example from Equation ( 12) and in various multi-peaked examples, conjecturing that such minima ought to approximate asymptotic solutions to the original problem in Equation (1).We then test these conjectures by investigating Equation (1) numerically.
which is a solution to Hence Using Equation ( 13) and rearranging gives Since πσ < √ 2L (see above Equation 12), this is a negative quadratic in ϵ.Furthermore, the maximum is where ϵ = p 2L .Now, ϵ ∈ 0, p 2L , so E[u * ] is an increasing function of ϵ on the interval 0, p 2L .Hence the minimum energy is where ϵ = 0.This analysis suggests that if a numerical solution to either Equation (1) or ( 5) results in a single peak at long times, we might expect that peak to be of a similar form to Equation ( 12) with ϵ = 0. We test this conjecture by solving Equation (1) numerically with initial conditions given by Equation ( 12) for various different values of ϵ ∈ 0, p 2L , fixing p = L = 1.For these simulations, we set D = 1, γ = 10, and so that σ = δ/ √ 3. Numerics reveal that the system does indeed tend towards a singlepeaked solution, where the width of the peak is approximately √ 2πσ and the solution is zero elsewhere (Figure 2).However, the asymptotic distribution is more flat-topped than the initial condition, owing to the fact that the initial condition arises from a moment closure approximation of K * u.This approximation reduces the analytic solution to a single Fourier mode, whereas the numerical solution could have arbitrarily many Fourier modes.
Finally note that, in the case ϵ = 0.3 (Figure 2c,f), a second peak emerges around x = ±1 (which are identified due to the periodic boundaries, recalling that L = 1).However, this decays by about t = 4.We will return to this phenomenon of decaying secondary peaks in the next section.

Twin peaks
In this section, we examine situations where there are two peaks.First, we look at situations where the peaks are the same height, then at cases where one peak is smaller than the other.

Peaks of identical height
Similar to the single-peak case, here we want to understand whether it is energetically favourable for a solution to have no mass outside the two peaks.More precisely, we examine the energy of the following solution to Equation (5), which is a critical point Fig. 3 Similar to the single peak case (Figure 2), when we start with two peaks of equal heights, surrounded by an area of constant density ϵ, that area becomes sucked-up into the peak.Panels (a) and (b) show this for x 0 = 0.5 and ϵ = 0.2, where both peaks remain.For x < 0.5, peaks merge, shown in Panel (c) for x 0 = 0.25.Panel (d) shows the time to merge as a function of x 0 .Parameters D, γ, and K are as in Figure 2. of otherwise. (20) Here, is half the (shortest) distance between the centres of the two peaks.As in the single-peak case, we can use Equation ( 7) to calculate A direct calculation using the definition of E[u] from Equation ( 8) leads to Since √ 2πσ < L, this is a negative quadratic in ϵ.The unique turning point is a maximum at ϵ = p 2L , so E[u * ] is an increasing function of ϵ on the interval 0, p 2L .Hence the minimum energy in the two-peaked case is where ϵ = 0, as with the onepeaked case.However, comparing the ϵ = 0 situation with one peak (Equation 18), against that with two peaks (Equation 22), we see that the single peak is a lowerenergy solution.This suggests that we might also see a merging of the two peaks, as well as the mass outside the peaks tending to zero.Indeed, in our numerical experiments, we saw a merging of peaks except in the special case where x 0 = 0.5, so that the initial peaks are evenly-spaced.Figure 3a,b shows an example where x 0 = 0.5 but ϵ > 0.Here two peaks remain but the the mass outside those two peaks is absorbed into the peaks over time.Figure 3c gives an example of peak merging for x 0 < 0.5 whilst Figure 3d shows how the time it takes for peaks to merge increases dramatically as x 0 increases towards x 0 = 0.5.Here, the time to merge is defined as the time at which the centre of the two initial peaks drops below 0.1.Whilst this is a rather arbitrary definition, other definitions lead to similar trends.
Notice here that the energy analysis does not give direct insight into why merging does not happen for x 0 = 0.5.Instead, we turn to physical intuition: the fact that peaks are evenly-spaced means that there is no 'preferred' direction for them to move in order to coalesce.Therefore they remain as two peaks.

Peaks of differing heights
In Section 4.1, we examined situations where there are two peaks with precisely equal height, finding that both peaks persisted indefinitely when they are evenly-spaced.However, we have already seen in Figure 1 that when peaks are of different heights, the smaller one can shrink over time, whereas the larger one grows.If this continues indefinitely, the smaller peak could decay completely and only one peak would remain, although it might take a long time for this to happen.
Here, we seek to explain this phenomenon using our energy approach, ascertaining whether we should always expect a smaller peak to end up decaying to zero, or whether there are situations where two peaks remain.To this end, we examine steady state solutions with the following functional form In this case, we can use Equation ( 7) to calculate We see immediately that, in order for c A and c B to be non-negative, we must have c A , c B ∈ 0, p √ 2πσ .A direct calculation using the definition of E[u] from Equation (8) leads to This is a negative quadratic in c B with critical point at c B = c A = p 2 √ 2πσ .Therefore the energy minima occur either when In other words, they occur when there is just one peak.Consequently, away from the critical point where c A = c B , we would expect the smaller peak to slowly decay to zero over time, leaving just one peak.Indeed, this is what we see in numerical solutions of Equation ( 5) (e.g. Figure 4a,b).However, the time it takes for the smaller peak to decay can be very large (Figure 4c).This is exacerbated by decreasing the diffusion constant, D (Figure 4d).Here, numerics hint that, as D → 0, the decay time may  19) with δ = 0.1.The value of c A is determined by Equation ( 24).
tend to infinity, meaning that the second peak may persist indefinitely if there is no diffusion to allow the smaller to seep into the larger.

Including population growth
So far, we have studied a system where the population size remains constant.This assumes that there are negligible births or deaths on the timescales that we are studying.Our focus has been on examining the difference between long transients and asymptotic solutions.However, in any real biological system, the effect of births and deaths will become non-negligible at some point in time.Therefore there is a limit to which transient solutions in these systems are biologically realistic: if the transients persist for too long, it will become necessary to account for the effect of births and deaths in any biologically meaningful model.We therefore examine the extent to which incorporating growth might enable a second peak to persist, by solving the following equation numerically with initial conditions given by Equation ( 23).
Depending upon the values of γ, D, K, and c B , we found that there is a critical value r = r c above which the second hump persists, and below which it decays.Figure 5a,b shows this in the case γ = 10, D = 1, K = 5, c B = 1, whereby r c ≈ 0.23. Figure 5c demonstrates how r c depends upon the aggregation strength γ: the greater the aggregation strength, the higher the required growth rate to enable a second peak to persist.

Discussion
Distinguishing between asymptotic solutions and long transients in numerical PDEs is a thorny issue, with perhaps no one-size-fits-all solution.Typically, researchers decide that a solution has reached an asymptotically-stable state when some measure (e.g. the change in L p norm for some p ∈ [1, ∞]) is below a small threshold value (see e.g.Burger et al (2014); Giunta et al (2022a); Schlichting and Seis (2022)).However, this means that if transient solutions are changing slower than this threshold value then they will be mistaken for asymptotically-stable solutions.Therefore it is valuable to have some analytic insight to guide the user as to whether the solution is (or is likely to be) a long transient or an asymptotically-stable solution.
Here, we have provided such a deductive technique for the aggregation-diffusion equation in Equation (1).Rather than studying this equation directly, we instead study an approximation given in Equation ( 5).This approximate formulation is simple enough to solve for steady state solutions.It also possesses an energy functional, which allows us to search for local minimum energy solutions amongst the steady state solutions, an approach employed successfully in a previous multi-species study (Giunta et al, 2022b).Our hypotheses are first that these local minimum energy solutions are stable solutions to Equation (5), whereas other steady states are not; and second that this categorisation carries over to the steady states of Equation (1).In the examples we tested, numerical experiments confirmed these hypotheses, with the sole exception of twin-peaked solutions where the peaks are of identical height and evenly-spaced.We therefore conclude that this method is a useful way for guiding users (i.e.those wanting to solving Equation 1 numerically) as to whether a solution they are observing is likely to be stable or not, whilst also recommending that they verify these calculations up with numerical experiments.
Regarding the examples we tested, we found two main results: first, that stable aggregations are likely to resemble compactly-supported solutions, rather than being non-zero everywhere; second, that multi-peaked solutions will always be transient unless either D = 0 or the peaks are precisely the same height and evenly-spaced.In addition to these central messages, further numerical investigations revealed that these twin-peaked transient solutions can be arbitrarily long-lived if the peaks are arbitrarily close to being evenly-spaced (Figure 3) and the heights of these peaks are arbitrarily similar (Figure 4).
That said, the consideration of very long transients in a model that operates on timescales where births and deaths are negligible is not terribly realistic, so we also examined the effect of adding a small amount of (logistic) growth.We found that arbitrarily small amounts of growth will not stop the smaller peak from decaying.However, there appears to be a critical growth rate, dependent upon the model parameters, below which the smaller peak will decay and above which it will grow (Figure 5).Therefore, if long transients appear when using Equation (1) to model biological aggregation, it is valuable to think about the effect of net reproductive rate in the system being modelled, and whether this is sufficient to arrest the decay of the smaller peak.
Whilst our principal equation of interest is Equation (1), it is worth noting that our approximate analytic techniques can also be applied to various other Equations.For example, the cell adhesion equations introduced in Armstrong et al (2006), have a very similar functional form that can usually be formally related to Equation (1) or modifications thereof (Painter et al, 2023).Chemotaxis equations are also somewhat similar to Equation ( 1), but here the non-local self-interaction is replaced with a diffusing chemical.The organisms interact with the chemical rather than directly with one another.It turns out that the resulting models are equivalent to a type of aggregationdiffusion equation with advection that is nonlocal in both space and time (Shi et al, 2021).This contrasts with Equation (1), which is nonlocal in space alone.However, similar patterns are observed in these systems, including long-transient multi-peaked solutions similar to those studied here (Potapov and Hillen, 2005).We also note that the moment closure we use in Equation (1) leads to a fourth-order PDE quite similar in nature to the Cahn-Hilliard equation (Novick-Cohen, 2008), for which there is a long history of studies on metastability (Bates and Xun, 1994;Reyna and Ward, 1995;Scholtes and Westdickenberg, 2018).
Finally, it is worth noting that the particular version of the aggregation-diffusion equation that we study involves linear diffusion.However, there is also interest in the nonlinear case, particularly where the diffusion is quadratic, replacing u xx with (u 2 ) xx = 2(uu x )x in Equation ( 1).An advantage of this formulation is that Equation (1) has the form u t = [u(D − γK * u) x ] x , making it amenable to a analysis without taking the limit D/γ → 0. This fact has been exploited, for example, by Ellefsen (2021); Carrillo et al (2018).However, here we have chosen to focus on linear diffusion is important to study as it often arises naturally from models of organism movement (Armstrong et al, 2006;Potts and Schlägel, 2020;Painter et al, 2023).Future work on the nonlinear case could reveal analytic insights about the effect of D vs. γ on asymptotic patterns, which we were only able to examine numerically in this study.

Fig. 2
Fig.2When the initial condition is a single peak surrounded by an area of constant density ϵ, that area becomes sucked-up into the peak.Panels (a) and (d) show this for ϵ = 0.1; (b) and (e) have ϵ = 0.2; (c) and (f) have ϵ = 0.3.In the latter case, a second peak emerges at x = ±1 but decays by around t ≈ 4, to leave a single-peaked final state.Panels (a-c) show the time-evolution of the system.Panels (d-f) show the initial conditions (blue curves) and final states (black).In all panels, D = 1, γ = 10, and K is a top-hat kernel (Equation19) with δ = 0.1.

Fig. 4
Fig. 4 Panel (a) shows a numerical solution of Equation (5) with initial condition given by Equation (23) with c B = 1.5.Panel (b) gives snapshots of the initial and final distributions.Notice that the smaller peak has decayed almost completely by t ≈ 15.Panel (c) is constructed from numerical solutions of Equation (5) with initial condition given by Equation (23) but with c B taking a variety of values, giving different start heights for the smaller peak (note that the start height is 2c B ). Panels (c) and (d) plot the time it takes for the smaller peak to decay to a maximum height of less than 0.1.This increases exponentially as a function of the start height, explaining the appearance of longtransient multi-peaked solutions to Equation (5) (Panel c).Conversely, the decay time decreases as D is increased, showing how diffusion can speed up decay of the smaller peak (Panel d).In panels (a-c), D = 1.In panel (d), c B = 1.In all panels, γ = 10 and K is a top-hat kernel (Equation19) with δ = 0.1.The value of c A is determined by Equation (24).

Fig. 5
Fig. 5 Effect of growth parameter.Panels (a) and (b) show the initial condition (blue) and solution at time t = 10 (black) where the parameters are γ = 10, D = 1, and K = 5.In Panel (a), r = 0.23 whereas Panel (b) has r = 0.24.This demonstrates a transition in long-term patterns, whereby the smaller peak decays for r ≤ 0.23 but grows for r ≥ 0.24.Panel (c) shows how this transition point, rc, decreases exponentially as the strength of attraction, γ, increases.