Physical interactions in non-ideal fluids promote Turing patterns

Turing’s mechanism is often invoked to explain periodic patterns in nature, although direct experimental support is scarce. Turing patterns form in reaction–diffusion systems when the activating species diffuse much slower than the inhibiting species, and the involved reactions are highly nonlinear. Such reactions can originate from cooperativity, whose physical interactions should also affect diffusion. We here take direct interactions into account and show that they strongly affect Turing patterns. We find that weak repulsion between the activator and inhibitor can substantially lower the required differential diffusivity and reaction nonlinearity. By contrast, strong interactions can induce phase separation, but the resulting length scale is still typically governed by the fundamental reaction–diffusion length scale. Taken together, our theory connects traditional Turing patterns with chemically active phase separation, thus describing a wider range of systems. Moreover, we demonstrate that even weak interactions affect patterns substantially, so they should be incorporated when modelling realistic systems.

Natural periodic patterns, ranging from nanocrystals [1], tissues [2], populations dynamics [3], to geophysical phenomena [4], are often explained by the seminal Turing mechanism [5][6][7][8].Turing patterns generally describe the spatial distribution of an activator and an inhibitor that diffuse in space.Patterns then form when the localized activator triggers production while the inhibitor suppresses production globally, often summarized as local activation, global inhibition [6,7].However, it is not * L.M. and C.L. contributed equally to this work.
clear whether Turing's mechanism can actually explain natural patterns [9,10] since inhibitors need to diffuse much faster than activators and the involved reactions need to be highly non-linear [11][12][13].Such non-linear reactions are often motivated by co-operative reactions, where multiple reactants lead to non-linearities [12].Cooperativity typically originates from physical interactions, which should also affect the diffusive motion of the species, but this is typically not taken into account.
Physical interactions are crucial for organizing biomolecules in cells [14], cells in tissues [15], and even organisms in groups [16].In particular, multivalent interactions can induce phase separation, where a dense droplet phase segregates spontaneously from a dilute surrounding phase.This is possible since the enthalpic gain from the interactions overcompensates the entropic loss of concentrating constituents [17,18].In simple passive systems, surface tension implies that large droplets grow to system size at the expense of small droplets [19].However, chemical reaction can suppress this Ostwald ripening [20] and thus control the size and arrangement of droplets [21][22][23].The predicted hexagonal arrangements [24] are very reminiscent of Turing patterns, although patterns are driven by phase separation in these systems.Taken together, although droplets regulated by chemical reactions share some properties with Turing patterns, it is unclear how the two models are related.
In this paper, we study a minimal system that is capable of forming droplets as well as Turing patterns.Effectively, we add physical interactions between activator and inhibitor to a standard reaction-diffusion system.This approach allows us to quantify how interactions affect pattern formation and it unveils a range of systems that interpolate between Turing's mechanism and patterns formed in active phase separating systems.In particular, we show how weak repulsive interactions stabilize patterns by inducing cross-diffusion, while strong interactions lead to phase separation, where coarsening is arrested by chemical reactions.

A. Interactions affect pattern formation
We start by considering a minimal pattern forming system comprised of two species: an activator A and an inhibitor I.The basic Turing model describes the dynamics of the respective fractions φ A (r, t) and φ I (r, t) as a function of the spatial position r and time t, for i = A, I. Here, the first term on the right hand side describes ideal diffusion with a diffusivity matrix D ij and the second term captures chemical reactions based on the Hill-Langmuir equation [25].For h ≥ 1, these reactions promote the production of A and I by activator A and suppress it by the inhibitor I, while both species exhibit linear degradation with rate k.This choice of chemical reactions allows us to independently control the typical fraction φ 0 of components A and I, the reaction rate k, and the reaction non-linearity h.We show in the Supporting Information that Eq. 1 exhibits a Turing instability for sufficiently large h and diffusivity ratio D I /D A , so the inhibitor spreads out while the activiator stays localized.
To include physical interactions between activator A and inhibitor I, we first consider the thermodynamics of an incompressible, isothermal fluid comprising the species A and I as well as an inert solvent S; see Fig. 1A.This system is still fully described by the volume fractions φ A and φ I , since the solvent occupies the remaining fraction φ S = 1 − φ A − φ I .The interactions of A and I in such a fluid can then be described by the Flory-Huggins free energy [26][27][28] where the integral is over the volume of the system, k B T is the relevant energy scale, and ν denotes a molecular volume, which we assume to be the same for all species.
The first three terms in the square bracket capture the translational entropies of all species, the fourth term describe the physical interaction between A and I, and the last term limits the width of interfaces between coexisting phases to roughly w in strongly interacting systems [28].The interactions between A and I are quantified by the Flory parameter χ: Positive χ denotes effective repulsion, which can originate from heterotypic repulsion or homotypic attraction, while negative χ leads to attraction of A and I. Equilibrium states, which minimize the free energy given by (2), can be inhomogeneous when interactions are sufficiently strong [29]: For strong attraction (large negative χ), a phase enriched in A and I will segregate from one enriched in the solvent, whereas strong repulsion (large positive χ) will lead to segregation of A from I with an equal amount of solvent in both phases.However, it is unclear how this equilibrium behavior is modified by the active reactions described by (1) and how Turing patterns are affected by weak interactions.
To model reaction-diffusion dynamics with interactions described by the Flory parameter χ, we replace the ideal diffusion term in Eq. 1 by a more general form which describes diffusion in non-ideal fluids.Linear non-equilibrium thermodynamics implies that diffusive fluxes are then proportional to gradients of the chemical potentials associated with the free energy given by Eq. 2, and the proportionality constants (known as Onsager coefficients or mobilities) determine the kinetic rate [30,31].Defining non-dimensional exchange chemical potentials we thus find where D i are the diffusivities of the species i = A, I, which are related to the mobilities D i φ i in this multicomponent system [32].We show in the Supporting Information that Eq. 4 reduces to Eq. 1, and thus describes ideal diffusion, if physical interactions are absent (χ = 0) and the wave length of patterns is large compared to w.Consequently, (4) describes a reaction-diffusion system encompassing non-ideal diffusion and containing normal Turing patterns as a limiting case.
To see how interactions affect patterns, we performed numerical simulations of Eqs.(3)-( 4) in a onedimensional system with periodic boundary conditions; see Methods.Fig. 1 demonstrates that without interactions (χ = 0), patterns with finite amplitudes only emerge if the reactions are sufficiently non-linear (large h) and the inhibitor diffuses sufficiently fast (D I D A ), as expected for Turing patterns [5].This trend persists for weak interactions, although the corresponding threshold values of h and D I /D A change.Apparently, repulsion between A and I promotes pattern formation (χ > 0), while attraction suppresses it (χ < 0).However, very strong attraction can again lead to large amplitudes (χ −9), independent of h and D I /D A .
The corresponding volume fraction profiles shown in Fig. 1D corroborate these observations: Without interactions (middle column), the system stays either homogeneous (pink parameter set) or forms normal Turing patterns (red parameter set) with a localized activator A and a fairly homogeneous inhibitor I.In contrast, strong attraction (left column) leads to co-localization of A and I, reminiscent of phase separation, albeit with a well defined pattern length scale.Similarly, strong repulsion (right column) implies segregation of A from I. Taken together, we thus showed that there is an interesting interplay between stereotypical Turing patterns and interactions promoting phase separation.

B. Weak interactions imply cross-diffusion
To understand how interactions affect pattern formation, we first analyze weak interactions (|χ| 5) by treating them as a perturbation to normal Turing patterns.Assuming the wave length of patterns is large compared to w, the generalized diffusion in Eq. 4 can be approximated by ideal diffusion to first order in χ; see Supporting Information.Consequently, the dynamics are described by Eq. 1 with the diffusivity matrix where ψ = φ 0 /(1−2φ 0 ).This analysis demonstrates that without interactions (χ = 0) in dilute system (φ 0 1, avoiding crowding effects), diffusion is dominated by the diagonal entries, resulting in stereotypical Turing patterns.In this case, we show analytically in the Supporting Information that patterns can only form when the diffusivity ratio D I /D A and the reaction non-linearity h are sufficiently large, consistent with the literature [12,13].
If A and I interact (χ = 0), Eq. 5 reveals that interactions directly affect cross-diffusion of A and I.For example, repulsive interactions (χ > 0) imply fluxes of A opposite to the gradient of I, thus favoring the segregation of the two species and enhancing patterns [33].To quantify this behavior, we analyze the covariance, cov(φ A , φ I ) = φ A φ I − φ A φ I , where the brackets denote spatial averages in the stationary state.Fig. 2 shows that the covariance generally decreases with more repulsive interactions (larger χ, consistent with enhanced cross-diffusion), increasing reaction non-linearity h, and diffusivity ratio D I /D A .The more detailed stability analysis presented in the Supporting Information demonstrates that repulsive interactions always promote pattern formation and lower the required reaction nonlinearity h and diffusivity ratio D I /D A ; see Fig. 3.In contrast, attractive interactions (χ < 0) generally stabilize the homogeneous system.However, this behavior only holds for moderate interactions χ since strong attraction (χ −9) also leads to large amplitudes; see Fig. 1.

C. Strong interactions invoke phase separation
Turing's mechanism cannot explain patterns that form when the activator A and inhibitor I attract each other strongly (χ < −9).Based on the strong correlations between A and I seen in Fig. 2, we hypothesize that strong attraction leads to associative phase separation of A and I from the solvent, while the reactions play a minor role.Indeed, we show in the Supporting Information that phase separation is possible in the absence of reactions when χ < χ − , where χ − = 8 arctanh(1 − 4φ 0 )/(4φ 0 − 1) marks the binodal point for a given average fraction φ 0 of A and I. Fig. 1 shows that the value χ − is very close to the onset of patterns, and that the resulting profiles are perfectly co-localized.Taken together, strong attraction between A and I leads to co-segregation of the two components from the solvent while the reactions barely affect the amplitude.
Strong repulsion between A and I should also lead to phase separation.In fact, for χ > χ + with χ + = 1/φ 0 , we predict that A can segregate from I spontaneously, even without reactions present; see Supporting Information.Figs.1-2 suggest that increasing the repulsion beyond this points indeed results in strong anti-correlation between A and I and a vanishing threshold for h and D I /D A .This numerical data indicates a continuous transition from patterns formed by reactions and diffusion (Turing patterns for weak interactions) to those formed by phase separation (strong interactions, χ < χ − or χ > χ + ).Taken together, we showed that patterns can form by reactions and by phase separation with an intricate interplay between them.

D. Reaction rate controls pattern length scale
We next ask what determines the length scale of the patterns.We show in the Supporting Information that is hardly affected by variations of the reaction nonlinearity h and diffusivity ratio D I /D A .In contrast, the interaction strength χ has a stronger influence: More repulsive interactions lead to patterns with shorter wave lengths, presumably because larger χ promote pattern formation.However, the strongest influence on the pattern length scale is the reaction rate k: Numerical simu-lations and the linear stability analysis presented in Fig. 4 indicated that k allows adjusting over several orders of magnitude with barely any changes in the pattern amplitude.
To understand how the reaction rate k affects the pattern length scale , we first focus on weak interactions.In this case, interactions mainly cause cross-diffusion (see Eq. 5 and Supporting Information), implying that the reaction diffusion lengths D A /k and D I /k are the only length scales in the equations.Consequently, length scales in the stationary state and in the initial instability must scale with k −1/2 for weak interactions, consistent with Fig. 4B.For strong interactions (χ < χ − or χ > χ + ), the system exhibits phase separation, implying that the initial instability is dominated by short patterns of length w while the stationary state patterns may exhibit much longer length scales due to coarsening [34].Fig. 4B shows that the linear stability analysis indeed predicts ∼ w in the region where we predict phase separation.For associative phase separation at strong attraction (χ < χ − ) these patterns remain stable and coarsening is suppressed; the variation in reflects the influence of χ on the interfacial width; see Supporting Information.In the contrasting case of strong repulsion (χ > χ + ), patterns coarsen to the reactiondiffusion length and thus scale with k −1/2 ; see Fig. 4B.This behavior is similar to the coarsening observed in active droplets, where the final length scale is also governed by the reaction-diffusion length [20].Note that the numerical data presented in Fig. 4B might not represent the full stationary state since domain sizes only grow logarithmically with time in these 1D systems [34].However, our analysis demonstrates that the length scale of the final pattern is generally governed by reaction-diffusion lengths, except when strong attraction between A and I leads to associative phase separation.

E. Results generalize to higher dimensions
So far, we have focused on pattern formation in one dimension for simplicity, but many natural patterns form in planar geometries.To see whether our results hold for this relevant case, we next perform a few selected simulations in two dimensions; see Fig. 5. Analogously to one dimension, we find Turing patterns for weak interactions (middle column of Fig. 5A) and strong interactions induce phase separation.In particular, strong attraction between A and I leads to co-localization (left column) whereas strong repulsion induces anti-correlated patterns (right column).Interestingly, in both cases of phase separation droplets form instead of stripe patterns, even though both phases occupy roughly half of the space.In such a case, normal phase separating systems exhibit stripe patterns [35], but the reaction-diffusion dynamics in our system apparently alter the picture.In any case, Fig. 5B shows that the pattern length scales we measured in 1D are very close to the ones measured in 2D simulations, implying that the results from the simple 1D system translate to the more complex 2D system.

III. DISCUSSION
We propose an extension to Turing patterns that take into account physical interactions that occur naturally.Weak repulsion between activator and inhibitor enhances patterns by inducing cross-diffusion, thus amplifying local activation and global inhibition.In contrast, strong interactions lead to phase separation, which can either be associative (A and I co-localize) or segregative (A separates from I).Both cases exhibit patterns for a much larger range of diffusivities and reaction non-linearities than normal Turing patterns, and the resulting length scales differ: In the segregative case of strong repulsion, patterns are governed by the reaction-diffusion length scale and thus grow larger for weaker reactions.In contrast, patterns in the associative case of strong attraction are arrested at the interfacial width.Taken together, we thus demonstrated that interactions can affect patterns substantially.The linear stability analysis presented in the Supporting Information demonstrates that these results do not depend on the specific choice of the reactions.Instead, interactions can generally lift restrictions on diffusivities and reaction non-linearities imposed by ordinary Turing patterns.Since physical interactions are virtually always present, many natural patterns can likely be explained by similar mechanisms.
Physical interactions in natural systems can stem from various sources and are virtually unavoidable in multicomponent systems.We need to investigate such systems in more detail, both in terms of physical interactions [30], chemical reaction networks [11,36], and conservation laws [37].For instance, Turing patterns can form when two species have equal diffusivity, while a third one is immobile [38], to produce effective differences in diffusivities.Explaining natural patterns in detail also requires incorporating growth [39], flows [40], noise, and delays [10].Moreover, natural patterns often form in complex geometries, including coupled layers [41] and curved surfaces [42], where the mechano-chemical coupling [43] can lead to dynamic patterns [44].The organization of biological cells is a particularly exciting example since biomolecules are known to interact and react [14].While this sometimes leads to spatial patterns explained by Turing's mechanism [7] other examples are akin to active droplets [21].Another possibility is patterns formed by self-propelled agents, which can exhibit motility-induced phase separation [45], and explain some population patterns successfully [16].In all these cases, physical interactions will affect patterns qualitatively and quantitatively, opening new perspectives on how natural patterns emerge.

A. Methods
We perform numerical simulations of Eq. 4 on an equidistantly discretized grid using second-order finitedifferences to approximate differential operators [46].We evaluate ∇µ on a staggered grid to ensure material conservation and use an explicit Euler scheme for the time evolution.

B. Acknowledgement
We thank Pierre Haas, Riccardo Rossetto, Noah Ziethen, and Yicheng Qiang for helpful discussions and critical reading of the manuscript.We gratefully acknowledge funding from the Max Planck Society and the Eu-ropean Union (ERC, EmulSim, 101044662).

Appendix A: Linear stability analysis of full model
We here present details of the linear stability analysis of the dynamical equations, given by Eq. 4 in the main text, for the volume fractions for the activator, φ A (r, t), and the inhibitor, φ I (r, t).The homogeneous stationary state is given by φ A = φ I = φ 0 .We linearize Eq. 4 in the main text around the uniform solution and determine the time evolution of perturbations in Fourier space, where the perturbations are characterized by the wave vector q.Their stability is then determined by the eigenvalues of the Jacobian matrix [8] J where we defined ψ = φ 0 /(1 − 2φ 0 ), which is the ratio of the fraction of species A and I to the fraction of the solvent.The two eigenvalues of J (q) are given by For an instability to occur, the real part of one of the eigenvalues must be positve.The real part of σ − (q) is always negative since tr J (q) = − q 2 (D A + D I ) 1 + ψ + φ 0 w 2 q 2 + 2k < 0. In contrast, the real part of σ + (q) is positive when det J (q) < 0. To homogeneous state is thus unstable if there is at least one value of q for which this is the case and the bifurcation lines in Fig. 1-4 in the main text can be determined numerically solving the condition min q det J (q) = 0 . (A3) Moreover, we can determine the length scale = 2π/q max associated with the most unstable mode, where the rate σ + (q max ) is maximal.Eq.A2 shows that the maximum of J (q) coincides with the minimum of det J (q), allowing us to determine using numerical minimization.Although the linear stability analysis only provides information at the initial stage when the volume fractions are close to φ 0 , the final length scale of the pattern will often be similar to the instability length scale; see Fig. 4 in the main text and Fig. S1.

Appendix B: Limit of weak interactions
To compare our extended model with physical interactions (Eq. 4 of the main text) to ordinary Turing systems (Eq. 1 of the main text), we show a formal reduction.

Expansion of full model for weak interactions
We thus consider the limit of weak perturbations, where the entropic contributions of the solvent (third term in the integrand in Eq. 2 of the main text) and the gradient terms (last term in the integrand) are negligible.In this case, we find that the diffusive terms in Eq. 4 reduce to ideal diffusion with the diffusivity matrix where ψ = φ 0 /(1 − 2φ 0 ).The dynamics in this limit are described by Eq. 1 in the main text with D as the diffusivity matrix.In particular, the stability of the homogeneous state φ A = φ I = φ 0 follows from the eigenvalues of the matrix which is the same as J (q) given in Eq.B2 except for the q 4 terms.The bifurcation lines can again be determined using Eq.A3.Fig. S2 shows that the resulting bifurcation lines are very similar to those determined for the full model, particularly for χ ≈ 0, thus validating the approximation.

Linear stability analysis of approximate model
We next determine bifurcation lines by solving Eq.A3.The determinant of J (q) is a quadratic polynomial in q 2 , det( J (q)) = a(q 2 ) 2 + bq 2 + c (B3) with where η = D I /D A is the diffusivity ratio.To find the minimum of det( J (q)), we can distinguish two cases based on the sign of a.For a < 0, the minimum of det( J (q)) is at q 2 = +∞ and always negative.The condition a < 0 reduces to χ < χ − spin or χ > χ + spin , where consistent with the analysis of the spinodal lines presented in Section C. In the opposing case of a > 0, the wave vector associated with the most unstable mode obeys d(det J (q))/d(q 2 ) = 0 and thus reads q 2 max = −b/(2a), or, If −b/(2a) < 0, we have det( J (q)) > det( J (0)), implying that the homogeneous state is stable since det( J (0)) = k 2 .However, if −b/(2a) > 0, the minimum value of det( J (q)) is min q det( J (q)) = det( J (q = q max )).In order to make the uniform state unstable, we thus require Consequently, the homogeneous state is unstable if where u = (η − 1)η(η(h − 1) + h + 1).Taken together, we get three bifurcation curves in the parameter space that we are interested in: (i). χ = χ − spin = −1/(φ 0 (1 − 2φ 0 )) for all h > 1 and η > 1.The uniform state is unstable for χ < χ − spin .
Fig. S2 summarizes these results and shows that the stability analysis of this limit of weak interactions (red lines) is very close to the full numerical stability analysis (black lines).Furthermore, Fig. S3 shows the critical value χ * as a function of h and η = D I /D A .For weak interactions (white region, χ * ≈ 0), both decreasing h and decreasing D I /D A leads to higher value of χ * , which agrees with the special cases shown in Figs. 1 and 3 in the main text.Another interesting observation is that the chemical reaction rate k does not appear in Eq.B9, implying the stability-instability transition is independent of the specific value of k, consistent with the results shown in Fig. 4 of the main text.

Analysis of most unstable wave length
We next use the approximate model in the limit of weak interactions to analyze the pattern length scales shown in Fig. 4 of the main text and Fig. S1.Here, = 2π/|q max | is associated with the most unstable mode given by Eq.B6.In particular, we find (1 − φ 0 χ) 2 + 1 + η (1 + 2ψ + φ 0 χ) 2 > 0 (B12) if η > (h + 1)/(h − 1).Consequently, the length scale also decreases as χ increases.All these variations of capture the trends observed in Fig. 4 of the main text and Fig. S1.In the limit η → ∞, we furthermore find implying the length scale converges to a constant value for a given k and h.In contrast, for given k and η, we have for large h, so increasing the reaction non-linearity can in principle result in arbitrarily short .
For segregative phase separation (χ > 0), the activator and the inhibitor repel each other.Assuming the amount of A and I are equal and the interactions are symmetric, the solvent fraction is identical in both phases, implying φ which now depend on φ and φ 0 .To determine the minimal value of the interaction where (spontaneous) phase separation can be expected, we minimize both expressions with respect to φ, which results in χ + spin = χ + bin = 1/φ 0 .
In the main text, we show binodal points for φ 0 = 0.2, resulting in χ − = χ − bin = −8.11and χ + = χ + bin = 5.Fig. S4 demonstrates that the width of equilibrium interfaces is of the order of the microscopic length w and increases with proximity to the binodal lines given by χ ± .The scaling exponent − 1 2 is consistent with the behavior close to the critical point [28].

FIG. 1 .
FIG. 1. Interactions affect pattern formation.(A) Schematic of chemical and physical interactions of activator A, inhibitor I, and the inert solvent.(B) Stationary state amplitudes of fraction φA as a function of the interaction strength χ and the reaction non-linearity h for diffusivity ratio DI /DA = 5. (C) Amplitude as a function of χ and DI /DA for h = 5. (D) Stationary patterns of φA (blue) and φI (orange) for the indicated parameters.(B, C) The homogeneous state is stable between the white lines, obtained from a linear stability analysis of Eq. 4, and the gray triangles mark critical interaction values χ− and χ+; see Supporting Information.(B-D) Model parameters are k = 0.1 DA/w 2 and φ0 = 0.2.Simulations ran for t = 10 5 w 2 /DA on a 1D grid of length 200 w with periodic boundary conditions.

FIG. 2 .
FIG. 2. Interactions affect correlation between activator and inhibitor.Covariance φAφI − φA φI as a function of interaction χ and reaction non-linearity h (left, DI /DA = 5) or diffusivity ratio DI /DA (right, h = 5).The homogeneous state is stable between the black lines, obtained from a linear stability analysis of Eq. 4; see Supporting Information.Model parameters are k = 0.1 DA/w 2 and φ0 = 0.2.Simulations ran for t = 10 5 w 2 /DA on a periodic 1D grid of length 200 w.

FIG. 3 .
FIG. 3.Repulsive interactions improve trade-off between differential diffusivity and reaction nonlinearity.Minimal interaction χ * to support patterns (Eq.S11 in the Supporting Information) as a function of diffusivity ratio DI /DA and reaction non-linearity h for φ0 = 0.2.Turing patterns without interactions (χ * = 0) form above the solid line.Conversely, phase separation is required to form patterns below the dotted line (χ * > χ+ = 5).

FIG. 4 . 1 0
FIG. 4. Reaction rate k determines pattern length scale.(A) Amplitude of activator φA as a function of interaction χ and reaction rate k.(B) Pattern length scale determined from the maximum of the factor of φA from numerical simulations (left) and from the fast growing mode in a linear stability analysis (right) as a function of χ and k (C) Stationary patterns of φA (blue) and φI (orange) for various parameters indicated in panels A and B. (B, C) The homogeneous state is stable between the white lines, obtained from a linear stability analysis of Eq. 4, and the gray triangles mark critical interaction values χ− and χ+; see Supporting Information.(A-C) Model parameters are h = 5, DI /DA = 10, φ0 = 0.2, and k0 = DA/w 2 .Simulations ran for t = 10 7 k −1 0 on a periodic 1D grid of length 2000 w.
FIG. S1.Pattern length scale determined from the maximum of the structure factor of φA determined from numerical simulations (left) and predicted by a linear stability analysis (right).The upper panels show as a function of interaction χ and diffusivity ratio DI /DA for reaction nonlinearity h = 5.The lower panels show as a function of χ and h for DI /DA = 10.Remaining model parameters are k = 0.1 DA/w 2 and φ0 = 0.2.Note that length scales do not vary much as a function of h and DI /DA, while χ has a significant influence.