Cross-diffusion induced patterns for a single-step enzymatic reaction

Several different enzymes display an apparent diffusion coefficient that increases with the concentration of their substrate. Moreover, their motion becomes directed in substrate gradients. Currently, there are several competing models for these transport dynamics. Here, we analyze whether the enzymatic reactions can generate a significant feedback from enzyme transport onto the substrate profile. We find that this feedback can generate spatial patterns in the enzyme distribution, with just a single-step catalytic reaction. However, patterns are formed only for a subclass of transport models. For such models, nonspecific repulsive interactions between the enzyme and the substrate cause the enzyme to accumulate in regions of low substrate concentration. Reactions then amplify local substrate fluctuations, causing enzymes to further accumulate where substrate is low. Experimental analysis of this pattern formation process could discriminate between different transport models.


Introduction
Experiments performed during the last decade found that at least eight different enzymes display a higher diffusion coefficient when the concentration of the corresponding substrate 1 arXiv:2006.02195v1 [physics.bio-ph] 3 Jun 2020 in solution is increased [1][2][3][4][5][6][7] . These increases are in the range of 24-80% relative to the diffusion coefficients without substrate. Most of these experiments relied on fluorescence correlation spectroscopy (FCS) measurements. Although artifacts introduced by this technique have been pointed out 8,9 , recent findings using other techniques 10,11 validated the phenomenon, which is often referred to as "enhanced diffusion". The underlying mechanism is still under debate 12 . Some experiments suggest that catalysis plays a key role 3,5,13,14 , while others indicate that enhanced diffusion persists when the substrate is replaced by an inhibitor 1,4,15 .
According to the latest experiments performed with the enzyme urease 11, 13 , it appears that both the binding and the catalysis step of the reaction scheme contribute to enhanced diffusion.
A related question is how enzymes behave in the presence of substrate gradients. The answer to this question appears to be complex. Some experiments suggest that enzymes drift downstream gradients of substrates, performing "antichemotaxis" 5, 13 . Others suggest that enzymes move upstream gradients of substrates, performing "chemotaxis" 6,7,15 . Antichemotaxis can be explained based on just the enhanced diffusion 5, 16 (enzymes accumulate in regions with low substrate concentration where they have a lower diffusion coefficient), chemotaxis cannot be generated by enhanced diffusion alone. However, cross-diffusion is a possible cause for enzyme chemotaxis 6 . Cross-diffusion describes the response of the enzyme to forces generated by gradients of substrate. Mathematically, it corresponds to an off-diagonal element in the diffusion matrix describing the combined motion of the enzyme and the substrate. It has been suggested that cross-diffusion can be due to specific interactions (ligand binding) between the enzyme and the substrate 6,15,17 or due to nonspecific interactions (e.g. steric, van der Waals) 18 . Specific interactions only lead to chemotaxis, while nonspecific interactions can cause the enzymes to move both up-or downstream the substrate gradient, depending on whether the interactions are attractive or repulsive, respectively. The model including nonspecific interactions 18 can be considered as a mathematical generalization of other existing models 5,6,17 .
While the existing models and experiments study how enzymes move in pre-imposed substrate gradients, they do not consider the feedback from the enzymatic reaction onto the substrate distribution. Here, we analyze the effects of this feedback starting from the most general transport model 18 . We show that spatial patterns can emerge in initially homogeneous systems if nonspecific interactions contribute to the accumulation of the enzyme in regions where concentration of substrate is low. Enzymes accumulating in these regions further deplete the substrate, causing the substrate gradient to become steeper, hence further increasing the accumulation of the enzyme. We obtain a set of conditions for the parameter range in which patterns form. We see that patterns can emerge only for repulsive nonspecific interactions, hence only for the model proposed in Ref. 18 , but not for the models proposed in Ref. 6,15 and Ref. 5 , suggesting that the analysis of pattern formation experiments can be used to discriminate between the different proposed models. Our findings imply that patterns can arise for a single-step enzymatic reaction even in the absence of autocatalytic activity or allosteric regulation. This is surprising given that the formation of conventional Turing patterns 19 with such simple reaction schemes requires at least a reaction network of three states, where forward and backward reactions are catalyzed by two different enzymes 20 .

Results
Our model system is depicted in Fig. 1. We consider a single-step enzymatic reaction in a narrow reaction chamber connected via a permeable membrane to a large substrate reservoir. We assume only substrate and product molecules can diffuse through the permeable membrane (exchange rates s , p , respectively), while enzymes are confined to the reaction chamber (since enzymes are typically larger than their substrates and products). The reservoir has a fixed concentration of substrate R and no products. For the reaction occurring in the bulk of this effectively 1D system we assume a Michaelis-Menten scheme (Fig. 1, Figure 1: System considered in this work. A single step enzymatic reaction takes place in a narrow reaction chamber which is coupled to a reservoir through a permeable membrane. The membrane allows for the exchange of only substrates and products with rate s , p respectively, but no exchange for the enzymes. Inset: Reaction scheme. The substrate is converted into a product by an enzyme ℰ. The reaction follows a Michaelis-Menten scheme where substrate can bind to the free enzyme with rate on , forming a complex, and unbind with rate off . The catalytic step of the reaction has a rate cat . We assume that the enzyme has a diffusion constant f , c depending on whether the enzyme is free or in its complexed form. Nonspecific pairwise interactions fs , cs can depend on the enzyme form. unbind with rate off . The catalytic step of the reaction has rate constant cat and catalysis is irreversible. This leads to a turnover rate per enzyme cat ( ) := cat /( M + ), where is the substrate concentration and M = ( off + cat )/ on . Following Ref. 18 we assume that the diffusion coefficient of the enzyme depends on whether the enzyme is free, f , or it is in its complexed form, c . This change in diffusion coefficient can be due to changes in either the hydrodynamic radius or the conformational fluctuations of the enzyme upon substrate binding 4,21,32 . Typically, enzymes become more tightly folded upon substrate binding, i.e. c > f , consistent with the experimentally observed trend. We also assume that short-range nonspecific interactions can either cause the enzyme to move towards or away from substrate, effectively generating a phoretic drift velocity that can also be interpreted as cross-diffusion 18 . We denote these interactions via pairwise potentials fs , cs depend-ing on whether the enzyme is free or in the complexed form, respectively, as depicted in Fig. 1

(inset).
Under the assumptions that (i) the enzyme is very dilute and (ii) the system is locally in chemical equilibrium (i.e. the timescales of diffusion and cross-diffusion are slower than the chemical reactions), one can derive an effective transport equation for the enzymes 18 , within the quasi-1D reaction chamber ( describes how enzymes respond to gradients of substrate due to the short-range nonspecific and hydrodynamic interactions, Here, c/f = Derjaguin length is a parameter capturing the effective short range interaction between the complex/free enzyme and the substrate. It is typically a few angstroms 24,25 , which is smaller than the Debye length (screening length) in typical buffer conditions (≈ 1 ) 12 . It is ex- 2 c/f is positive (negative) when the interaction is attractive (repulsive) 18 . The derivation of 2 c/f is similar to that of the second Virial coefficient for a real gas 26 , but it also includes hydrodynamic corrections and is computed by assuming that the size of the enzyme is much larger than the interaction length. The sign of 2 c/f determines the sign of c/f and ultimately the sign of xd ( , ). For attractive interactions xd ( , ) < 0, i.e. the enzyme drifts towards higher concentrations of substrate (chemotaxis). The enzyme performs antichemotaxis for repulsive interactions.
Note that the effect of nonspecific interactions can also be written as a phoretic drift 18 by swapping the in Eq. (1) with in the definition (2), with a drift velocity directly proportional to the substrate gradient ph ( , In the regime where the enzyme is dilute and Eq. (1) is valid, the reaction chamber of where ( , ) denotes the product concentration. While the product dynamics does not influence the substrate and enzyme equations, we include it as spatial read-out of the reaction.
Having the reaction chamber coupled to a reservoir avoids product accumulation and substrate depletion, generating a nonzero homogeneous steady-state, with concentrations h , h , and h for enzyme, substrate, and product, respectively. For simplicity, we express h and h as functions of h , Since ( h ) is a monotonic increasing function of h , it is possible to write h and h in terms of h , which can be directly tuned in experiments via the total enzyme concentration (see SI  Table S1 for the full param-  Table S1. In Fig. 2A we see that the homogeneous solution is unstable and patterns form, for a value of 2 depending on the parameters, the system will spontaneously form patterns .
To characterize the instability of the homogeneous steady state solution, the spatial frequency of the linear perturbation and the perturbation growth rate. If > 0 the perturbation grows with time and patterns form. A positive can be found provided that (see SI) with = s R . This result holds in the strong depletion regime ( h ≪ R ), where the effect of the reaction on the substrate dominates over the outflow to the reservoir. In this regime, the expressions are simpler and it is easier to pinpoint the driving mechanism behind the instability observed in Fig. 2. We refer the reader to the SI for the full analysis. The argument of the square root in the inequality (7) is positive if relation (6) is fulfilled. In the SI we show that the instability is a Type II instability 27 , meaning that = 0 at = 0. This is natural as the total amount of enzymes in our system is conserved and homogeneously increasing or decreasing perturbations, i.e. perturbations at = 0, would correspond to changes in the total enzyme amount. By looking at the inequalities (6), (7), we can see that repulsive interactions, i.e. c/f < 0, are needed to have instabilities. Diffusion tends to homogenize the concentration profiles and contributes with a positive term to the left hand side of inequality (6). In the case where the enzyme-substrate interaction is attractive, i.e. f > 0, also the second term on the left hand side of inequality (6) is positive. Hence the interaction between the enzyme and the substrate would need to change sign upon substrate binding to have an instability ( c < 0). There can be cases in which the change in nonspecific interactions is less abrupt and both f < 0 and c < 0. Even in the case for which c = f = , i.e. there is no change in interaction upon substrate binding, it is possible to have an unstable homogeneous solution for where we rewrote relation (6) (3) is unstable and patterns arise, the ratio˜/˜can be as high as 10 9 . Above the curve, the homogeneous solution is stable and the ratio˜/˜= 1. The white cross corresponds to the simulation shown in Fig. 2A. For the full parameter list we refer the reader to the SI.
It is interesting to note how relations (6)-(8) depend on the substrate concentration. One could ask, given certain nonspecific interactions between the enzyme and the substrate, at which substrate concentration * h should we begin to observe instabilities? From relations (6)- Above the instability lines the homogeneous solution is stable. Below the line the system is unstable and patterns similar to the one shown in Fig. 2A Figure 4: Positive feedback mechanism behind the pattern forming process. (i) We start with a substrate gradient and a homogeneous enzyme profile; (ii) the enhanced diffusive current e D and the cross-diffusive one e xd cause the enzyme to accumulate where substrate is low; (iii) reaction causes the substrate gradient to get steeper; (iv) this leads to a further increase of e D and e xd causing further accumulation of the enzyme. This process repeats itself determining the patterns.
What is the physical mechanism underlying the instability given by the inequalities Both e ( ) and ( ) have a Michaelis-Menten dependence on for all the models proposed so far for the enzyme motion. They differ only in the prefactors and a non-zero offset for e ( ). Inequality (9) is never fulfilled for such e ( ) and ( ) and consequently any initial perturbation of the concentrations is smoothed out by diffusion.
However, the inequality (9) will apply to any model with a spatially dependent diffusion coefficient that is coupled to another diffusing and reacting species. Systems of this type are The patterns observed in Fig. 2 are not generated via the common short-range activation and long-range inhibition mechanism 33 , as neither the enzyme nor the substrate have autocatalytic activity. It is also not a motility induced phase separation (MIPS) mechanism 34 , which relies on the slowing down of active particles in regions of high particle concentrations. Here a positive feedback mechanism between particle accumulation and reaction leads to the pattern formation. Note that for the system of equations (3) patterns for a simple enzymatic reaction. We found that enhanced diffusion would need to be more sensitive to perturbations in substrate concentrations than the reaction, see Eq. (9).
Although this is not the case for enzymes, we believe that relation (9)

Linear Instability Analysis
In this section, we perform the linear instability analysis of the PDE system (3) of the main text to derive the condition for which patterns form. As a reminder, the starting model for our analysis (Eq (3) of the main text) is:

Homogeneous steady state solution
The homogeneous steady state solution ⃗ v h = ( h , h , h ) of (10) is given by the following expressions: where h , h , and h are the homogeneous concentrations of enzyme, substrate, and product respectively. In this work for simplicity we consider h and h as functions of h but it is also possible to write h and h in terms of h , which is a quantity directly tunable in the experiments: which and In

Instability of the well-mixed system
We first study the instability of the solution for the well-mixed system. In the well-mixed system, the spatial derivatives and any spatial dependence of the concentrations are neglected.
The system (10) becomes: The fixed-point of the system (17) is still given by ⃗ v h . To perform a linear instability analysis, we perturb the system around its fixed point. Any perturbation around ⃗ v h can be written where ⃗ v( ) = ( ( ), ( ), ( )). The amount of enzymes in the system is fixed and does not change over time because ( ) = 0, hence any perturbation ( ) just shifts the total enzyme amount. The linearized system of equation takes the form: We rewrite the linear system of equations (18) in matrix form where The (10) and we linearize the system by assuming that perturbations are small and find that: where ′ e ( h ), ′ ( h ) are the derivatives of e ( ), ( ) with respect to at = h We write the linear system of equations (22) in Fourier space to partially diagonalize the equations. The result is that the dynamics of each Fourier mode is independent from other modes and is governed by the following equation: where and It is not surprising that J(0) is the same as the Jacobian of the well-mixed system J 0 , (21).
At = 0 we are neglecting all the effects due to diffusion and cross-diffusion, cf. Eq. (21).
Moreover perturbations with = 0 corresponds to homogeneous shifts in the concentrations, which are the same as considered for the well-mixed system.
Using the definition of xd , e and Eq. (11) we can write the following By rewriting the above expression as ( − 1 )( − 2 ) = 0, we see that the summation of the eigenvalues is negative, therefore the smaller eigenvalue 2 is always negative. The larger eigenvalue 1 is positive if and only if the product of the two eigenvalues is negative.
The condition for having negative product can be written as a condition for the wave vector In order to have a positive *2 , we should have that: Equations (30) and (31) are the instability conditions for the system of equations (10) (Eq. (3) of the main text). In the regime h ≪ R , with = s R , we get: which are identical to the relations (4) and (5) Figure 6: The largest eigenvalue of the Jacobian 1 versus the wave vector . The horizontal black line corresponds to 1 = 0. By increasing the initial homogeneous concentration of substrate, h , from below to above critical concentration, * h , the eigenvalues become positive. The shape of the eigenvalue curves indicates that the instability is of type II. The parameters for the graph are = 300 , f = 10 2 / , c = 13 2 / , s = 100 2 / , M = 1 , value, * h , which is related indirectly to the initial amount of enzymes in the system. The form of curves in Fig.6 correspond to a type II instability 27 . A type II instability is typical of systems with conserved quantities. In our case the total amount of enzymes is conserved, this implies that 1 = 0 at = 0, otherwise we would have homogeneous change in h over time, corresponding to changes in the total amount of enzymes.

No short range interactions
In this section we consider xd ( , ) = 0. The Jacobian (27) then takes the form: 3 = − p 2 − p is one of the eigenvalues and it is negative. The other two eigenvalues are obtained by solving the following equation: By rewriting the above expression as ( − 1 )( − 2 ) = 0, we see that the sum of the two eigenvalues is negative, as ′ ( h ) > 0. The largest eigenvalue is positive if the product of the two eigenvalues is also negative. Thus the term with zeroth order of should be negative.
Together with the identity Eq. (11), i.e. cat h = s ( R − h )/ ( h ), we find the following condition: In the regime h ≪ R , with = s R , we get: which is the same as the inequality (7) of the main text.
By using the definition of e ( ) used in 18 e = f + ( c − f ) ( ), we find that the inequality (40) takes the form: Because R > h and ′ ( h ) > 0, the left hand side of the inequality above is negative and cannot be larger than the always positive right hand side. Therefore enhanced diffusion alone cannot drive instabilities. The cross-diffusion generated by repulsive interactions is key in this respect. This also holds in the simpler case of h ≪ R , with = s R .
Similarly it is possible to show that we cannot have instabilities for other enhanced diffusion definitions 5, 15 .

Parameters of the simulations
The numerical simulations of the system of equations 10 have been carried out by using the COMSOL Multiphysics v5.3.
We used the parameters listed in Table 1 for the simulations shown in Fig. 2  concentrations h , h , h were perturbed by adding white Gaussian noise with variance h /100, h /100, h /100 respectively. We considered periodic boundary conditions and we used the "Time dependent" solver of COMSOL with a relative tolerance of 10 −6 and an element size of 0.01 , i.e. 10 4 lattice points.
For the results shown in Fig 3 we again used the "Time dependent" solver of COMSOL.
We considered a grid of 20 × 20 values for 2 c = 2 f = 2 in the interval [−10, 3]Å 2 (evenly spaced), and h in [10 1 , 10 5 ] (evenly spaced on a log-scale). For the COMSOL specific parameters, we used a relative tolerance of 10 −6 and a maximum element size of 0.02 , except when two peaks were observed as a final result of the simulation. In these cases, we repeated the simulations with a finer grid (element size of 0.01 ) and again we observed a single peak for the concentrations as the final result of the simulations. All the other parameters were the same as of Table 1 and the boudary conditions were periodic.