Pattern formation is an emerging phenomenon observed across various disciplines, including biology, chemistry, physics, and engineering. Examples range from the stripes on a zebra to the spirals of galaxies [1]. Understanding the mechanisms and principles behind pattern formation has been a longstanding scientific challenge, as it plays a fundamental role in the organization and functionality of natural systems. In biology, pattern formation underpins key processes such as embryonic development, tissue morphogenesis, and cellular differentiation, which are essential for the proper structure and function of living organisms. Mathematical models play a crucial role in understanding pattern formation, offering a powerful tool to unravel its intricate dynamics, uncover fundamental principles, and predict complex behaviors [2][3][4]. During past decades, significant contributions have advanced our understanding of pattern formation, including the molecular underpinnings [5], modeling dispersal in biological systems [6], investigating dynamic pattern formation in cellular networks [7], and exploring multiscale models of developmental systems [8].
Reaction-diffusion equations, rooted in the seminal work of Alan Turing [3], have been one of the best-known mathematical models used to study pattern formation [9][10][11][12]. Arising from Turing instability, reaction-diffusion models can generate a wide variety of spatial patterns. Recent advances in the mathematical study of pattern formation include the integration of computational techniques to compute multiple patterns [13][14][15][16][17][18], as well as the application of machine learning methods to analyze complex pattern formation from data and discover novel patterns [16,19].
One of the most well-known reaction-diffusion systems for pattern formation is the Gray-Scott model, which can generate a diverse range of spatial patterns [20][21][22][23]. Originally proposed by Gray and Scott [24,25], this model consists of two reacting chemical species in a system of two ordinary differential equations. Later, mechanical effects, i.e., diffusion, was incorporated into the model, leading to reactiondiffusion equations capable of generating complex patterns [20]. The Gray-Scott model contributes to our understanding of pattern formation observed in various natural and artificial systems [9] has been widely studied during the past decades [23,26,27]. These insights not only deepen our understanding of nonlinear dynamics and self-organization but also have practical implications in fields such as chemistry, materials science, and biology [28]. The Gray-Scott model is a widely used prototype for studying pattern formation, providing a critical link between theoretical modeling and experimental observations across various disciplines [29]. Similar models for pattern formation include the Belousov-Zhabotinsky reaction [30], the Schnakenberg model [31], the Gierer-Meinhardt model [32], and the Brusselator model [33], which also exhibit rich dynamics and are commonly used in mathematical biology and related fields.
The classical Gray-Scott model is built by coupling linear diffusion equations with the reaction kinetics of two irreversible reactions, described by the law of mass action. More precisely, the Gray-Scott model considers the following two irreversible reactions
where 𝑉 is an activator, 𝑈 is a substrate, and 𝑃 is an inert product. To maintain the system out of equilibrium, both 𝑈 and 𝑉 are removed by the feed process [20]. Consequently, The reaction-diffusion equation of 𝑢 and 𝑣 are given by [20]
subject to certain boundary and initial conditions. Here, 𝑢 and 𝑣 denote the concentration of 𝑈 and 𝑉 , the dimensionless reaction rate for the first reaction is set to be 1, 𝑘 is the dimensionless rate constant of the second reaction, 𝑓 represents the dimensionless feed rate, and 𝐷 𝑢 and 𝐷 𝑣 are diffusion coefficients. A surprising variety of spatiotemporal patterns can emerge from the reaction-diffusion Eq. ( 2) under specific parameter values and initial conditions [20].
From a modeling perspective, the classical Gray-Scott model is not thermodynamically consistent, meaning it may not satisfy the first and second laws of thermodynamics. In biological systems, where energy input, dissipation, and conversion govern self-organization, thermodynamic consistency is crucial for developing mechanistically accurate and predictive models that faithfully capture the underlying physical principles. As a result, the thermodynamic basis for pattern formation in the Gray-Scott model is unclear. For instance, it remains uncertain whether these patterns represent transient phenomena or nonequilibrium steady states, what is the energetic cost of the complex pattern in the models.
In a recent work [34], a variational, reversible Gray-Scott model is developed based on an energetic variational approach (EnVarA) [35,36]. The model revises the classical Gray-Scott model into a thermodynamically consistent form by including all reverse reactions in the model. Additionally, a virtual species is introduced to account for the influx and efflux of 𝑢 and 𝑣, effectively transforming the open system into a subsystem of a larger, closed system. It is proved in [34] that in the short-time dynamic process, when 𝜖, the small constant related to reaction rates in the reverse part, goes to zero, the solution to the variational Gray-Scott models will converge to that of the classical Gray-Scott models. Additionally, numerical simulations in a recent study [37] demonstrate pattern formation can occur in the variational model.
The variational Gray-Scott model not only restores thermodynamic consistency of the classical model but also provides a more realistic framework for studying biological and chemical pattern formation. Many natural processes, such as biochemical reactions and cellular transport mechanisms, involve reversible dynamics that adhere to thermodynamic laws. By capturing these aspects, the variational model allows for a deeper understanding of pattern formation in systems where classical models fall short. Moreover, the techniques and insights developed through this study can be extended to other reactiondiffusion systems, offering a roadmap for analyzing complex patterns in fields ranging from developmental biology to materials science.
For the classical Gray-Scott model, many nontrivial steady states on the domain 𝛺 = (0, 1) with no-flux boundary conditions have been computed in [13]. As a first step in exploring pattern formation in the thermodynamically consistent variational Gray-Scott model for various values of 𝜖, we use these steady states as initial conditions, expecting they may serve as suitable starting points for observing potential pattern formation in the variational model. Numerical experiments indicate that for relatively large values of 𝜖, all initial conditions quickly converge to a uniform steady state. However, for sufficiently small 𝜖, the variational model exhibits rich transient dynamics: stationary patterns may persist for long periods, and both oscillatory and traveling-wave-like behaviors can emerge. We also analyze the stability of two uniform steady states in the variational Gray-Scott model for fixed 𝜖 to provide some theoretical insights to the simulation results. Although both uniform states are stable in a certain sense, pattern formation can only occur if the initial condition does not converge to the boundary steady state, in which (𝑢, 𝑣) = (1, 0).
The rest of the paper is organized as follows. In Section 2, we derive the variational Gray-Scott model using the energetic variational approach. Section 3 is devoted to analyzing the energy stability of two uniform steady states. In Section 4, we investigate pattern formation in one dimension for different values of 𝜖, using the steady states from [13] as initial conditions. Finally, in Section 5, we examine how the persistence time of patterns depends on the parameter 𝜖. A conclusion remark is given in Section 6.
In this section, we briefly review the derivation of the thermodynamically consistent variational Gray-Scott model, proposed in [34], by an energetic variational approach [35,36].
Motivated by non-equilibrium thermodynamics, particularly the celebrated works of Lord Rayleigh [38] and Onsager [39,40], the EnVarA has been a powerful tool of building thermodynamically consistent models in physics, chemical engineering, and biology [35,41]. The key idea of EnVarA is to describe an isothermal and mechanically isolated system by its energy and the rate of energy dissipation over time, along with kinematic (transport) assumptions on the employed variables. More precisely, according to the first and second laws of thermodynamics [35,42], an isothermal and closed system possesses an energy-dissipation law
Here, 𝐸 total is the total energy, which is the sum of the Helmholtz free energy and the kinetic energy ; ▵ (𝑡) ≥ 0 stands for the rate of energy dissipation, which equals to the rate of entropy production in this case. Once these quantities are specified, a thermodynamically consistent model can be derived by combining the Least Action Principle (LAP) and the Maximum Dissipation Principle (MDP) [35]. More specifically, for the energy part, one can employ the LAP, taking variation of the action functional (𝒙) = ∫ 𝑇 0 ( - ) d𝑡 with respect to 𝒙 (the trajectory in Lagrangian coordinates) [35,43], to derive the conservative force, i.e., 𝛿 = ∫ 𝑇 0 ∫ 𝛺 (force iner -force conv ) ⋅ 𝛿𝒙 d𝒙d𝑡. For the dissipation part, one can apply the MDP, taking the variation of the Onsager dissipation functional with respect to the ''rate'' 𝒙 𝑡 , to derive the dissipative force, i.e., 𝛿 = ∫ 𝛺 force diss ⋅𝛿𝒙 𝑡 d𝒙, where the dissipation functional = 1 2 ▵ in the linear response regime [39]. Consequently, the force balance condition results in
which is the dynamics of the system. In the case that = 0, 𝛿 𝛿𝒙 = -𝛿 𝛿𝒙 , then the dynamics can be written as
which is a generalized gradient flow. For these systems, the free energy determines the equilibrium of the system and the rate of energy dissipation determines the dynamics. The EnVarA is originally developed for mechanical systems, and 𝒙 should be understood as the flow map. Recent work [41] extends this variational principle to reaction kinetics by computing the variation with respect to the reaction trajectory 𝑅, analogous to the flow map, and its time derivative 𝜕 𝑡 𝑅, representing the reaction rate, to derive reaction kinetics. One key advantage of using EnVarA to model complex chemo-mechanical systems in biology is its ability to provide a unified framework for both mechanics and chemistry. All multiscale energetic couplings and competing effects are naturally incorporated through the choice of the energy-dissipation law.
Since we are only interested in the concentration of the substrate 𝑈 and the activator 𝑉 , we can combine the removal process 𝑉 and the process of generating 𝑃 , an inert product, as one process for mathematical simplicity. Therefore, the chemical reactions described in the classical Gray-Scott model (2) can be represented as follows:
In the reaction network (5), all chemical reactions are irreversible. In addition, it involves the birth and death of 𝑈 , making the overall system an open system.
To formulate a thermodynamically consistent variational Gray-Scott model, we consider the following reversible chemical reactions:
where 𝜖 𝑖 are small parameters, 𝑌 is a virtual species added to model the birth-death process of 𝑈 . Clearly, if 𝜖 1 → 0 and 𝜖 2 → 0, the first two reactions revert to those in (5). The introduction of the virtual species 𝑌 , whose concentration is large (of the order 𝑂(1∕𝜖 3 )), allows us to treat the open system, situated in a sustained environment with influx and efflux, as a subsystem in a larger, closed ''universe'' [44]. A similar approach is used in [45]. In general, the small parameters 𝜖 𝑖 are different and may not be in the same order. In particular, 𝜖 3 not only represents the reaction rate of the reverse reaction but also describes the concentration scale of the virtual species 𝑌 . In the current study, we assume 𝜖 1 = 𝜖 2 = 𝜖 3 = 𝜖. We will explore more general cases in future work. We denote the concentrations of species 𝑃 and 𝑌 by 𝑝 and 𝑦∕𝜖 respectively. To derive the variational Gray-Scott model, we first notice that the concentrations satisfy the following kinematics
where 𝑅 𝑖 is the reaction trajectory of each reaction in (6), and 𝒖 𝑢 and 𝒖 𝑣 are effective velocity induced by the diffusion process. The reaction trajectories 𝑅 𝑖 (𝑡) account for the number of forward reaction has happened for 𝑖th reaction by time 𝑡 and may take negative values [41]. Here, we disregard the diffusion effects on 𝑃 and 𝑌 , as the slow diffusion of P and Y minimally impacts the dynamics of U and V, which are our main focus. The boundary condition for all species is the non-flux boundary condition, given by
which ensures the boundary term vanishes in deriving the force balance equation by the EnVarA. It is important to note that
for the kinematics (7) along with the non-flux boundary condition (8).
The conservation property ( 9) plays an important role in studying the steady state of the variational Gray-Scott model. Following the general framework of modeling reaction-diffusion systems [41], the variational Gray-Scott model can be formulated based on the energy-dissipation law
Here, (𝑢, 𝑣, 𝑝, 𝑦) is the free energy of the system, ▵ mech and ▵ cheme are the rate of energy-dissipation due the mechanical (diffusion) and chemical (reaction) parts respectively. The free energy (𝑢, 𝑣, 𝑝, 𝑦) is taken as
Here, ỹ = 𝑦 𝜖 represents the concentration of 𝑌 , 𝜎 𝑢 , 𝜎 𝑣 , 𝜎 𝑝 , 𝜎 ỹ denote internal energy of each species, which together determine the equilibrium of the system. Let (𝑢 𝑠 + , 𝑣 𝑠 + , 𝑝 𝑠 + , ỹ𝑠 + ) ∈ R 4 + be a positive equilibrium of the chemical reaction system (6), then 𝜎 𝑢 , 𝜎 𝑣 , 𝜎 𝑝 , 𝜎 ỹ satisfies
Since
we have
We can solve for 𝜎 𝑖 (𝑖 = 𝑢, 𝑣, 𝑝, ỹ) in terms of 𝑘, 𝑓 and 𝜖. Notice in Eq. ( 14), that there are four variables but only three equations, resulting in an overparameterized case. The analysis for the differential internal energy 𝜎 𝑖 follows a similar pattern. In this paper, we take
Next, we impose the rate of energy dissipation of the system, given by
and
)
To derive the dynamics from the energy-dissipation law, we apply the EnVarA to the mechanical and chemical parts respectively, which leads to 𝒖 𝛼 (𝛼 = 𝑢, 𝑣) and 𝑅 𝑖 (𝑖 = 1, 2, 3) such that the energy-dissipation law (10) holds. In order to get the equation of 𝒖 𝑢 , we first define the flow map 𝒙 𝑢 (𝐗, 𝑡) that is associated with 𝒖 𝑢 by the ODE
If we consider only the diffusion of 𝑢, 𝑢(𝒙, 𝑡) is determined by the flow map 𝒙 𝑢 (𝐗, 𝑡), making as a functional of 𝒙 𝑢 (𝐗, 𝑡). Applying the LAP and MDP with respect to 𝒙 𝑢 and 𝒖 𝑢 leads to the force balance equation (see [34] for the detailed calculations)
where 𝜇 𝑢 is known as the chemical potential for the species 𝑢. Similarly, we can derive that
For the chemical part, we need to perform the EnVarA with respect to 𝑅 𝑖 and 𝜕 𝑡 𝑅 𝑖 , which leads to ln
Recall that 𝜎 𝑢 , 𝜎 𝑣 , 𝜎 𝑝 , 𝜎 ỹ satisfies ( 14), we can show 𝜕 𝑡 𝑅 𝑖 satisfies the law of mass action. Indeed, from (20), we know that
The calculations for 𝜕 𝑡 𝑅 2 and 𝜕 𝑡 𝑅 3 are similar.
Combining the force balance equations ( 18)-( 20) with the kinematics (7), we end up with the following reversible variational Gray-Scott model
Remark 2.1. Although the Gray-Scott model is not directly associated with any specific biological system, the same approach can be applied to restore thermodynamic consistency in more realistic models of pattern formation in various biological systems, as these models are typically of the reaction-diffusion type [12]. As another representative example, consider the celebrated Schnakenberg model [31,46], given by
which models the sequence of reactions
where 𝑢 and 𝑣 denote the concentrations of species 𝑋 and 𝑌 , respectively, and 𝑎 and 𝑏 are the concentrations of 𝐴 and 𝐵, assumed to be constant. A thermodynamically consistent version of the Schnakenberg model can be constructed by restoring the reversible parts of the latter two reactions with small reaction rate 𝜖 and formulating the full system in terms of 𝑎, 𝑏, 𝑢, and 𝑣. The resulting reaction-diffusion system can be derived from the energy-dissipation law d d𝑡 ∫ 𝑢(ln 𝑢 -1) + 𝑢𝜎 𝑢 + 𝑣(ln 𝑣 -1) + 𝑣𝜎 𝑣 + 𝑎(ln 𝑎 -1) + 𝑎𝜎 𝑎 + 𝑏(ln 𝑏 -1)
)
where 𝑅 𝑖 are reaction trajectories for the three involved reactions, and 𝒖 𝛼 (𝛼 = 𝑢, 𝑣) are average velocities due to the diffusion of 𝑢 and 𝑣. The internal energies can be taken as
In this subsection, we show the first two equations in the variational Gray-Scott model (22) can be reduced to the classical Gray-Scott model (2) when 𝜖 → 0.
Assume the initial concentrations of U, V, P, and Y as 𝑢 0 , 𝑣 0 , 𝑝 0 , and 𝑦 0 ∕𝜖, where 𝑦 0 = 𝑓 . When 𝜖 is small, we notice that 𝑦(𝑡) ≈ 𝑦 0 = 𝑓 is a nearly constant function since 𝑦 𝑡 ≈ 0. Thus, as 𝜖 approaches zero, the first two equations in the variational Gray-Scott model formally converge to the classical Gray-Scott models, as we drop all terms with 𝜖 and replace 𝑦 by 𝑓 in the 𝑢, 𝑣 equation. The above argument can be made more rigorously [34]. Since both the equations of 𝑝 and 𝑦 are linear, we have
By plugging Eq. ( 24) into Eq. ( 22), we rewrite the equations of 𝑢 and 𝑣 as follows: ), where 𝑀 is an order one constant with respect to 𝜖, then the following limits hold:
and 𝑦 0 𝑒 -𝜖𝑡 → 𝑦 0 as 𝜖 → 0. Hence, formally, the variational Gray-Scott model ( 25) converges to the classical Gray-Scott model
when 𝜖 → 0. We emphasize that it might be difficult to prove the assumption max{‖𝑢‖ 𝐿 ∞ , ‖𝑣‖ 𝐿 ∞ } ≤ 𝑀. So, the above limit is formal. However, the numerical simulations in the next section show this assumption holds for all tested initial conditions. It is worth emphasizing that the formal convergence does not imply that the performance of the variational Gray-Scott models and the classical model is identical, as the system may exhibit singularity with respect to 𝜖. In the following sections, we will study the dynamical behavior of the variational Gray-Scott model for various 𝜖.
Unlike the classical Gray-Scott model (when 𝜖 = 0), which has only one uniform steady state (𝑢, 𝑣) = (1, 0), for a given 𝜖 > 0, the variational system admits two uniform steady states. More precisely, for a given initial conditions (𝑢 0 , 𝑣 0 , 𝑝 0 , 𝑦 0 ) and the domain 𝛺 = (0, 1), we define the total mass of all species 𝐶 as
where 𝑦 0 = 𝑓 . Assume that (𝑢 𝑠 , 𝑣 𝑠 , 𝑝 𝑠 , 𝑦 𝑠 ) is spatially homogeneous steady-states of the variational Gray-Scott model (22), then (𝑢 𝑠 , 𝑣 𝑠 , 𝑝 𝑠 , 𝑦 𝑠 ) satisfies
subject to the constraint:
Here, |𝛺| = 1 is the size of the domain. Solving (28) with the constraint (29), we obtain two spatially homogeneous steady-states: one is a boundary steady-state given by
and the other is an interior steady-state expressed as
where
Here, the term interior means the steady state is in the interior of the stoichiometric compatibility class [47], defined by
In contrast, the term boundary refers to a steady state that lies on the boundary of this compatibility class.
Notice that 𝜆 = (𝜖 -1 ) and 𝐶 = (𝜖 -1 ), we can estimate the two steady states as follows:
The concentration of 𝑌 is (𝜖 -1 ) in the boundary steady state and is (1) in the interior steady state. Since the initial condition of 𝑌 is (𝜖 -1 ), it can be expected that a significant time is needed if the system would like to reach the interior steady state for small 𝜖.
Remark 3.1. Note that when 𝜖 → 0, we have
) .
Hence, the boundary steady state in the variational Gray-Scott model corresponds to the uniform steady state of the classical Gray-Scott model, while the interior steady state is not related to any steady state of the classical Gray-Scott model.
In this section, we analyze the energy stability of two uniform steady states with fixed 𝜖. Although the variational Gray-Scott model comprises four species -𝑈 , 𝑉 , 𝑃 , and 𝑌 , the conservation law constrains the system's degrees of freedom to three, corresponding to three reaction trajectories 𝑅 𝑖 . Specifically, based on Eqs. (7), the perturbation in the stability analysis satisfies
where 𝛼 𝑖 ∈ R for 𝑖 = 1, 2, 3. Thus we define the perturbation manifold as follows:
where, 𝒔 𝑖 correspond to the perturbation along the reaction trajectory 𝑅 𝑖 . We use 𝐸 denote the free energy (11) without spatial integration throughout this section.
For the interior steady state 𝒖 2 , we can prove that it is a local minimizer of the free energy (11) based on the following proposition: Proposition 3.1. For any 𝜖 > 0, the interior steady state
is a local minimizer of the energy 𝐸 on the manifold 𝒖 0 + 1 , where 𝒖 0 = (𝑢 0 , 𝑣 0 , 𝑝 0 , 𝑦 0 ∕𝜖) is the initial condition.
Proof. Without loss of generality, we perturb 𝒖 2 along the direction of 𝒔 2 of the three-dimensional manifold 1 (a similar analysis applies to other directions) and obtain:
By the definition of 𝜎 𝑝 and 𝜎 𝑣 , we conclude that 𝛿 = 0 is the unique solution of d𝐸[𝒖 2 +𝛿𝒔 1 ] d𝛿 = 0, indicating the critical point. Furthermore, by computing the second variation, we have:
Since 𝑣 𝑠 2 , 𝑝 𝑠 2 > 0 for any 𝜖 > 0,
> 0, indicating that this steady state is a local minimizer in the 𝒔 1 direction, similarly for other directions. □ Next, we examine the energy landscape near the interior steady state and its evolution as 𝜖 approaches zero. We plot the energy landscape around the steady state in three directions on 1 shown in Fig. 1 for 𝜖 = 0.01 and 𝜖 = 0.0001 respectively. While the behavior of 𝒔 1 remains stable, for 𝒔 2 and 𝒔 3 , we observe a notable change as 𝜖 diminishes: the stable region contracts significantly. Essentially, the inflection point approaches the steady state in this direction.
The reason for this phenomenon is that 𝑢 𝑠 2 will be close to 0 when 𝜖 is small, as 𝑢 𝑠 2 is at the order of (𝜖). Computing the second variation in the 𝒔 1 direction, we have
Notice that
> 0 requires 𝛿 < 𝑢 𝑠 2 . Consequently, as 𝜖 decreases, the stability region around the local minimizer shrinks.
In summary, we have shown that the interior steady state remains stable for all 𝜖 > 0. However, as 𝜖 decreases, the stable region also becomes smaller. Nonetheless, it is still challenging for the system to escape from the steady state, as it would require 𝛿 > 𝑢 𝑠 2 , which would result in 𝑢 < 0, an impossible scenario. Therefore, even for very small 𝜖, if the initial data is not far away from the steady state, the system will still converge to this interior steady state. However, when 𝜖 is small, ensuring that the initial conditions are not too far from the steady state requires 𝑝 to be of order 𝜖 -1 , meaning species 𝑃 must have a significant presence.
For the boundary steady state, the stability analysis differs from that of the interior steady state. First, we observe that this is not a critical point of the free energy. Recall 𝒔 1 = (-1, 1, 0, 0) ∈ 1 , then:
This expression tends to infinity as 𝛿 → 0 + , indicating it is not a critical point of the free energy. Indeed, we can compute the local minimizer of 𝐸 is the three-dimensional perturbation space 𝒖 1 + 1 , which satisfies:
This system has a unique nonzero solution such that
Hence, 𝒖 2 is the unique minimizer of the free energy 𝐸 in the space 𝒖 1 + 1 . The reason that 𝒖 1 is a steady state of the system is that the mobilities for 𝑅 1 and 𝑅 2 , given by 𝜖𝑣 3 and 𝜖𝑝, vanish, due to the disappearance of 𝑉 and 𝑃 . Consequently, the first and second reactions cease to exist:
We can show that the boundary steady state is a local minimizer of the free energy 𝐸 in a one-dimensional perturbation space 2 ∶= span{(-1, 0, 0, 𝜖)}, i.e., a perturbation along the direction of the third reaction trajectory: Proposition 3.2. For any 𝜖 > 0, the boundary steady state
is a local minimizer of the free energy on the manifold 𝒖 0 + 2 , where 𝐮 0 = (𝑢 0 , 𝑣 0 , 𝑝 0 , 𝑦 0 ∕𝜖) is the initial condition.
Proof. The proof follows a similar structure to that of Proposition 3.1. □
Although the boundary steady state is only stable along the direction (-1, 0, 0, 𝜖) and not stable in the (-1, 1, 0, 0) and (0, 1, -1, 0) directions. It still behaves as a nearly local minimizer when 𝜖 is small and (1) amounts of U, V, and P, and (𝜖 -1 ) of Y present initially. More precisely, since the first and second reactions account for only a fraction of 𝜖 in the system, it is reasonable to consider the perturbation along the direction 𝒔 = 𝒔 3 + 𝜖𝒔 1 + 𝜖𝒔 2 = (-1, 0, 0, 𝜖) + (-𝜖, 0, 𝜖, 0), direct calculation reveals that:
Additionally, numerical results further support this assertion, as shown in Fig. 2. The boundary steady state is a local minimizer in the onedimensional space for small 𝜖. This is primarily because most reactions in this scenario are governed by the third reaction, which is stable for the boundary steady state. Consequently, for certain initial conditions, the system will converge to the boundary steady state when 𝜖 is small. We term it as virtual stability. The reason for labeling this as ''virtual'' is that the third reaction is artificial and virtual, as it does not exist in the classical Gray-Scott models defined by Eq. ( 2). Moreover, when 𝜖 is not very small, it can be noticed that the boundary steady state is not a local minimizer, even within the restricted one-dimensional setting. As a result, a small perturbation around it leads to convergence towards the interior steady state.
Overall, we observe that the interior steady state is stable, while the boundary steady state is considered virtually stable. The virtual stability arises when V and P disappear quickly in the system. In such cases, the initial terms of P and V occupy only a fraction of 𝜖 in the entire system, while the term Y dominates almost entirely. This scenario satisfies the conditions for virtual stability, leading the system to converge to the boundary steady state.
In this section, we explore the pattern formulation in the variational Gray-Scott model for different 𝜖 in one dimension.
For the classical Gray-Scott model (where 𝜖 = 0), stationary spatial patterns correspond to non-uniform steady states of the system [13]. Additionally, traveling patterns exist, which correspond to traveling wave solutions [27]. In [13], the authors compute steady states of the classical Gray-Scott model in the domain 𝛺 = (0, 1) subject to the non-flux boundary condition under the standard finite difference discretization. The parameter values used are 𝐷 𝑢 = 5 × 10 -4 , 𝐷 𝑣 = 2.5 × 10 -4 , 𝑓 = 𝑦 0 = 0.04, and 𝑘 = 0.065. It shows that, under these parameter values and the non-flux boundary condition, in addition to the trivial uniform steady state (𝑢, 𝑣) = (0, 1), the classical Gray-Scott model has 6 non-uniform linearly stable steady states and 16 linearly unstable nonuniform steady states, shown in Fig. 3. These steady states are computed using homotopy continuation techniques. Their linear stability is determined by analyzing the eigenvalues of the Jacobian matrix at the steady states. Specifically, the spectrum is computed numerically by solving the characteristic equation for the Jacobian. If any eigenvalue has a positive real part, the steady state is unstable, while if all eigenvalues have negative real parts, the steady state is stable.
To explore pattern formation in the variational Gray-Scott model for different 𝜖, we use these 23 steady states of the classical Gray-Scott model as the initial conditions for (𝑢, 𝑣). The initial conditions of 𝑝 and 𝑦 are taken as
We impose the same no-flux boundary conditions as in [13]. We then examine the evolution of these initial conditions for various values of 𝜖, with the expectation that they serve as suitable starting points for capturing pattern formation in the variational model. These initial conditions can also be interpreted as spatially heterogeneous perturbations of the two uniform steady states analyzed in the previous section. We adopt a semi-implicit method, which treats the reaction part explicitly and the diffusion part implicitly, for the temporal discretization. The spatial grid size is ℎ = 1∕256. Since the initial concentration of species 𝑌 is 𝑦 0 ∕𝜖, for smaller 𝜖, more 𝑌 exist and the total mass 𝐶, defined in (27), is also larger. We consider 𝜖 = 10 -2 , 10 -4 , and 10 -6 , corresponding to a relatively large 𝜖, an intermediate-sized 𝜖, and a significantly small 𝜖, respectively. visualized through 𝑢(𝑥, 𝑡) for 𝑥 ∈ (0, 1) and 𝑡 ∈ (0, 40). The temporal step size is set to 𝛥𝑡 = 0.5. Numerical results show that decreasing the temporal step size yields quantitatively similar solutions. The initial conditions in Fig. 4(a) correspond to linearly stable steady states of the classical Gray-Scott model, while those in Fig. 4(b) correspond to linearly unstable steady states.
The simulation result shows that, for relatively large 𝜖, all initial conditions, including the one with (𝑢, 𝑣) = (1, 0), converge to the uniform interior steady state 𝒖 2 , as the concentration of 𝑈 will decrease to around 0. As a result, any non-uniform patterns in the limiting system are destroyed in a relatively short time for a relatively large 𝜖. This result is consistent with the theoretical analysis presented in the previous section, as the interior steady state 𝒖 2 is a minimizer of the free energy and the boundary steady state is not virtually stable for a relatively large 𝜖. Unlike the case of 𝜖 = 10 -2 , the initial condition with (𝑢, 𝑣) = (1, 0) can converge to the boundary steady state. Additionally, there are four initial conditions that converge to the boundary steady state 𝒖 1 quickly. The result is consistent with the theoretical analysis that the boundary steady state is virtually stable for small 𝜖. Other initial conditions tend to converge to the interior steady state 𝒖 2 as the concentration of 𝑈 decreases. However, the convergence is very slow, and non-uniform patterns can persist for a long time. Non-stationary patterns are still observable at 𝑡 = 1000. For further illustration, Fig. 6 More interestingly, for 𝜖 = 10 -4 , we observe oscillating solutions during the time evolution. In the original paper on the Gray-Scott model [24], Gray and Scott demonstrate that the system exhibits chemical oscillations even without diffusion.
To illustrate the oscillation in the variational Gray-Scott model, we examine the solution with the initial condition corresponding to the linearly unstable solution 4 (fourth image in Fig. 5(b)) in detail. Fig. 7(a)-(b) show the evolution of 𝑢(0.5, 𝑡) and 𝑣(0.5, 𝑡), as well as the total mass of 𝑢 and 𝑣, respectively. The plots show the damped oscillations in the concentration of 𝑢 and 𝑣, which are similar to the phenomenon reported in Fig. 4 in [24] for the irreversible Gray-Scott model without diffusion. After the oscillation, the solution behavior is similar to the solution with the initial condition corresponding to the stable solution 3 (fourth image in Fig. 5(a)).
3. 𝜖 = 10 -6 : significantly small 𝜖 Next, we consider 𝜖 = 10 -6 , a significantly small 𝜖. Fig. 8 shows the simulation results for 𝜖 = 10 -6 , visualized by 𝑢(𝑥, 𝑡) for 𝑥 ∈ (0, 1) and 𝑡 ∈ (0, 1000). Again, the initial conditions in Fig. 8(a) correspond to linearly stable steady states of the classical Gray-Scott model, while those in Fig. 8(b) correspond to linearly unstable steady states.
Like the case with 𝜖 = 10 -4 , due to the virtual stability, the initial condition with (𝑢, 𝑣) = (1, 0) converges to the boundary steady state quickly. Additionally, more initial conditions will converge to the boundary steady state quickly, compared with the case of 𝜖 = 10 -4 .
Interestingly, for the initial conditions corresponding to six nonuniform, linearly stable steady states in the classical Gray-Scott model, the profiles of 𝑢 and 𝑣 remain unchanged throughout the evolution. In other words, these linearly stable stationary patterns in the classical Gray-Scott model can be stabilized as transient states in the variational model for a very long time when 𝜖 is significantly small. We can view these solutions as quasi-steady states or quasi-stable patterns as the (𝑢, 𝑣)-component of the solution is unchanged. Fig. 9 shows the concentrations of 𝑈 , 𝑉 , 𝑃 , and 𝑌 at 𝑡 = 0 and 𝑡 = 1000 for the initial condition associated with the non-uniform linearly stable solution 5. It can be observed that although the concentrations of 𝑈 and 𝑉 remain unchanged, the concentration of 𝑃 increases while 𝑌 decreases by the same amount. The effective dynamics of the whole system is to transform 𝑌 to 𝑃 . From the numerical experiments, one can expect the system will reach the interior steady 𝑢 2 at the end as the concentration of 𝑃 increases and 𝑌 decreases. However, since the initial concentration of 𝑌 is much larger than that of 𝑃 , a significant time is needed to reach the steady state. Consequently, the spatial pattern can be stabilized for a long time if the concentration of 𝑌 is large, i.e., 𝑦(𝑥) ≈ 𝑓 . Additionally, there are four initial conditions, corresponding to linearly unstable steady-state 5, 6, 9, and 10, which will initially evolve towards a quasisteady state and will remain unchanged in the (𝑢, 𝑣)-components for will determine the dynamics either transform 𝑌 to 𝑃 or 𝑃 to 𝑌 . In the current study, since the concentration of 𝑃 is taken as 1, which is much smaller than that of 𝑌 , the system will converge to the boundary steady state quickly if the essential dynamics is 𝑃 to 𝑌 . We can observe the pattern formation for a significantly long time if the essential dynamics is 𝑌 to 𝑃 . We will investigate the effects of the initial concentration of 𝑃 in future work.
Furthermore, the simulation results are consistent with the energy stability analysis presented in previous sections. When 𝜖 is large, all initial conditions converge to the interior steady state 𝒖 2 , which minimizes the free energy. However, when 𝜖 is sufficiently small, the gradient flow dynamics favor the boundary steady state 𝒖 1 for certain unstable initial configurations, due to the virtual stability of the boundary state in this regime. In this sense, both trivial steady states become effectively ''stable''. For some unstable steady states, the system tends to evolve towards the boundary steady state. Particularly, when the initial condition is close to the boundary steady state but far from the interior one -particularly due to the large value of 𝑌 -the system rapidly converges to the boundary state. In contrast, convergence to the interior steady state is much slower, leading to the persistence of non-uniform patterns.
In this section, we study the pattern persistence time in the variational Gray-Scott model with respect to 𝜖.
We consider two typical initial conditions, which correspond to the non-uniform linearly stable solution 6 (the last solution in Fig. 3(a)) and the linearly unstable solution 3 (the third solution in Fig. 3(b)). The first one will always converge to the interior steady state according to the numerical simulation. The second one will converge to the interior steady state for large 𝜖, but converge to the boundary steady state for small 𝜖. We define the pattern time as
with 𝑎 being significantly small since the vanishing of the pattern indicates that both 𝑢 and 𝑣 are close to constant functions. We pick 𝑎 = 0.05 throughout this section. Fig. 12 shows the relationship between the pattern persistence time and 𝜖 for these two initial conditions.
For the linearly stable solution 6, the pattern persistence time increases as 𝜖 decreases, and is of the order (𝜖 -1 ). For the linearly unstable solution 3, it converges to the interior state for relatively large 𝜖, and the pattern persistence time also increases with decreasing 𝜖, following the order (𝜖 -1 ). However, the initial condition converges to the boundary steady state for large 𝜖, and the initial pattern will be destroyed faster for smaller 𝜖. The result is consistent with the previous simulation and analysis.
We plot the evolution of free energy for two initial conditions with different values of 𝜖. For 𝜖 = 10 -2 , the free energy plot clearly shows that both initial conditions converge to equilibrium around 𝑡 = 400. The evolution of free energy for 𝜖 = 10 -4 follows a similar pattern, but it takes significantly longer for both solutions to reach the trivial steady state, as shown in Fig. 12. However, for 𝜖 = 10 -6 , the linearly unstable steady state in the classical Gray-Scott model quickly converges to the boundary steady state, although the free energy of this state remains relatively large. For the linearly stable state in the classical Gray-Scott model, while the (𝑢, 𝑣)-components of the solution remain unchanged, as shown in Fig. 13, the system's free energy continuously decreases over time. The free energy plot also suggests that the linearly stable To summarize, in order to have pattern formation in the variational Gray-Scott model, one needs to have some particular initial condition such that the variational model converges to the interior steady state, in which all 𝑌 will convert to 𝑃 . The pattern is maintained if the concentration of 𝑌 stays large, analogous to the continuous feed of 𝑈 in the classical Gray-Scott model. The conclusion is similar to that in a recent paper [50] on a slightly different thermodynamically consistent three-species reaction-diffusion model, in which the authors show that for a finite system, a specific Turing pattern exists only within a finite range of total molecule number, and the presence of the third species stabilizes the Turing pattern of the two species.
In this paper, we study the pattern formation of a thermodynamically consistent variational Gray-Scott model, derived by an energetic variational approach, in one dimension, using non-uniform steady states of the classical model computed in [13] as initial conditions. The variational Gray-Scott model includes a virtual term 𝑌 and reversible reactions to the classical Gray-Scott model, transforming the system into a thermodynamically consistent closed system. The classical Gray-Scott model can be viewed as a subsystem of the variational Gray-Scott model when the reverse reaction rate 𝜖 tends to zero. By decreasing 𝜖, we observed that stationary patterns in the classical Gray-Scott model can appear as transient states in the variational model when 𝜖 is significantly small. Additionally, the variational model admits oscillating and traveling-wave-like solutions for small 𝜖. The results show the capability of the variational model in capturing pattern formation.
We also analyze the energy stability of two uniform steady states, an interior steady state and a boundary steady state, in the variational Gray-Scott model. Although the interior steady state is always stable, the stability region becomes significantly smaller as 𝜖 decreases. In the meantime, the boundary steady state is virtually stable, i.e., is stable with respect to the third reaction 𝑈 ← ← ← ← ← ← ← ← ← ← ← ← ← ← ← ← ← ⇀ ↽ ← ← ← ← ← ← ← ← ← ← ← ← ← ← ← ← ← 𝑌 . For certain initial conditions, since the concentration of 𝑌 (at order (𝜖 -1 )) is much larger than other species (at the order of (1)) for small 𝜖, the gradient flow dynamics will drive the system to the boundary steady state. In order to observe pattern formation, one needs a special initial condition such that the dynamics will not converge to the boundary steady state.
The variational Gray-Scott model offers a new mathematical framework for understanding pattern formation from a thermodynamic perspective. The numerical simulation and theoretical analysis suggest that pattern formation is maintained by the presence of species 𝑌 , i.e., the continuous input of 𝑈 in the classical Gray-Scott model. The initial condition determines the effective direction of the reaction network, and the pattern formation can only occur if the network continuously generates the inert product 𝑃 . Furthermore, for sufficiently small 𝜖, the system's free energy decreases over time due to the conversion of 𝑌 to 𝑃 , indicating that these patterns act as transition states in a larger, closed system. This also implies that maintaining the pattern requires a continuous energy input from the environment. The current study represents a first step towards understanding pattern formation in biological systems from an energetic perspective. Several open questions remain for the variational Gray-Scott model, including: (1) Investigating the effect of varying the scale of small parameter 𝜖 𝑖 in the reaction network (6); (2) Analyzing the role of 𝑃 in pattern formation; (3) Examining more general initial conditions. We will study these open questions in future work. CRediT authorship contribution statement Wenrui Hao: Writing -review & editing, Supervision, Methodology, Funding acquisition, Conceptualization. Chun Liu: Writingreview & editing, Supervision, Methodology, Conceptualization. Yiwei Wang: Writing -review & editing, Writing -original draft, Software, Methodology, Investigation, Conceptualization. Yahong Yang: Writing -review & editing, Writing -original draft, Methodology, Investigation, Formal analysis, Conceptualization.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.