Entropy Production in Reaction–Diffusion Systems Confined in Narrow Channels

TL;DR

This paper studies how the shape of narrow channels affects entropy production in chemical reaction systems.

Contribution

The study introduces an effective diffusion equation to analyze entropy production in confined reaction–diffusion systems.

Findings

Entropy density changes in value but retains its qualitative behavior in narrow channels.

Channel geometry influences structure formation in the Gray–Scott system.

Abstract

This work analyzes the effect of wall geometry when a reaction–diffusion system is confined to a narrow channel. In particular, we study the entropy production density in the reversible Gray–Scott system. Using an effective diffusion equation that considers modifications by the channel characteristics, we find that the entropy density changes its value but not its qualitative behavior, which helps explore the structure-formation space.

Figures

Full text

1. Introduction

Many chemical and biological systems exhibit ordered structures that emerge spontaneously and remain stable over time. These structures arise and persist due to processes involving mass transport and chemical reactions under far-from-equilibrium conditions. For instance, Turing’s diffusion-driven instability mechanism describes how small fluctuations in reaction–diffusion systems can be amplified, forming spatio-temporal patterns [1,2]. These patterns arise from the interaction between chemical species, one considered inhibitory and the other activating. Diffusion induces the system’s uniform steady state, which is stable in the absence of diffusion, to become unstable, promoting the emergence of patterns. Transport originates in zones with higher inhibitor concentrations that block the activator, whereas low inhibitor concentrations do not, allowing the activator to grow [2]. This phenomenon is counterintuitive since diffusion is generally associated with a decrease in concentration gradients over time. The reaction–diffusion systems that have been used to elucidate pattern formation typically involve non-linear reactions. They often focus only on forward reactions, assuming their rates significantly surpass those of reverse reactions that have been used to elucidate these phenomena [3]. Consequently, for a complete understanding of these processes’ thermodynamics further study is needed. In classical thermodynamics, the evolution of systems within the linear regime is predictable. However, far-from-equilibrium systems may evolve into organized structures known as dissipative structures, as revealed by some analyses of irreversible thermodynamics [4,5,6].

Numerical studies have revealed some particular connections between entropy production and the emergence of patterns. Entropy production analysis in cycles of chemical oscillations or self-organizing biological systems has become possible in some reversible systems [7]. For instance, in chemical oscillations like the Belousov–Zhabotinsky reaction or the Oregonator model, entropy production averaged over a cycle may exceed or fall below the corresponding unstable equilibrium state, depending on the bifurcation parameter [8]. These processes are relevant, for instance, in self-healing surfaces involving pattern formation [9] or in enzymatic reactions of photosynthetic models that exhibit standard oscillatory and damped dynamics [10]. However, establishing thermodynamic rules for studying the evolution of non-equilibrium patterns remains challenging.

The Gray–Scott model is a powerful tool capable of replicating various Turing-like patterns. This model finds practical applications, accurately describing the autocatalytic reactions observed in chemical and enzymatic reactions in biochemistry and biology [11]. The original Gray–Scott model corresponds to two irreversible reactions for U and V, with P being an inert product [12]. This system models the chemical reaction $[eqn]$ , which consumes U and produces V. Therefore, the amount of both substances must be controlled to maintain the reaction. This is done by adding at the feed rate f and removing at the kill rate k. The reduction and subsequent removal of V can also be described by another chemical reaction $[eqn]$ . For this reaction, P is an inert product, meaning it does not react, and the parameter k controls the rate of the second reaction. The Gray–Scott system is defined by two equations that describe the behavior of two reacting substances:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ are the diffusion coefficients for each chemical species. In order to remove V faster than the growth of U, the rate k is added to f and multiplied by V, since it is assumed that the disappearance of V also depends on its concentration.

The reversible Gray–Scott system has been proposed to study entropy production [3], featuring a cubic autocatalytic reaction coupled with a first-order reaction. Escher, Ross [13], and Yoshida [14] also previously studied it in the context of continuous-flow well-stirred tank reactors (CSTR), where far-from-equilibrium situations can also be achieved. Its phase diagram presents a stationary-stable equilibrium point, saddle-node, and Hopf bifurcations. Heterogeneous structures typically arise near these bifurcation curves [15,16].

Entropy production in self-replicating pulses, initiated by perturbing the equilibrium point, has been investigated [3,17]. While entropy production increases with the size of each pulse, its rate gradually decreases. Thus, the entropy production value can serve as a thermodynamic metric for evaluating pattern stability. Additionally, total entropy change, incorporating outward entropy flux using a relative chemical potential, has been examined within this model [18,19]. Specific trends in entropy production evolution have been observed during the formation of heterogeneous structures. In contrast to uniform steady states, during this process a compromise occurs between maximization and minimization of entropy production, reflecting the influence of several interacting mechanisms [20].

Incorporating a thermal component into the system has revealed transitions between different morphologies of Turing patterns regulated by temperature [21]. These transitions occur with lower entropy production compared to other parameters, such as kinetic constants and reactant fluxes [22]. Formalisms such as the Boltzmann equation [23] and extended irreversible thermodynamics [24] have also drawn further insights into entropy production and entropy flows in reaction–diffusion systems.

On the other hand, the study of transport processes and reactions within confined domains holds significant interest due to its wide range of applications, including carbon nanotubes [25], zeolites [26], biosensors [27], biological membranes [28], separation processes [29,30], and microfluidics, among others [31]. Confinement plays a crucial role in influencing the transport properties by restricting the system’s accessible configuration space [32,33]. The Fick–Jacobs approximation for diffusion in narrow channels simplifies this process by assuming a fast equilibration of the transverse degrees of freedom. It focuses on the longitudinal motion by projecting along the channel axis, resulting in an effective one-dimensional equation with a position-dependent diffusion coefficient [33]. Numerous studies have proposed different functions of this diffusion coefficient under diverse conditions, such as in symmetrical channels [34,35,36,37], asymmetrical channels [38,39,40], channels with curved midlines [41,42,43], or on curved surfaces [44,45,46]. Additionally, confinement effects on chemical reactions [47] and pattern formation in reaction–diffusion systems [48], including diffusive predator–prey models [49] and long-range effects [50], were also explored. Significantly, the geometry of the channel walls introduces modifications to the Turing instability conditions, the dispersion relation, and the range of unstable modes [48].

This work studies how confinement in a narrow channel changes the entropy production density in a reaction–diffusion system. By introducing entropic flow into the Fick–Jacobs approximation, as in [51], it is clear that wall geometry affects entropy production. In the following section, we review the definitions of entropy production for reaction–diffusion systems, focusing on the reversible Gray–Scott model. Subsequently, we explore the effects of geometry on a confined reaction–diffusion system within the local Fick–Jacobs approximation and, particularly, within the mean-field effective framework. We conclude by presenting the effective entropy production calculated with both approaches and the results concerning entropy production density by diffusive flux.

2. Entropy Production in 2D Reaction–Diffusion Systems

Entropy production comprises two distinct components in reaction–diffusion systems, as discussed in [3,5,6,7,8]. The first component arises from chemical reactions’ entropy produced per unit area (or length). The second component is the rate of entropy production due to diffusion. This total entropy production is calculated by integrating the entropy production density across the entire spatial domain:

$$[eqn]$$

The contribution of chemical reactions is

$$[eqn]$$

where $[eqn]$ are the forward and backward reaction rates of the j-th reaction. Far from equilibrium, backward rates can be considered small, and the system’s affinity tends to infinity, leading to linear entropy production. In this scenario, chemical oscillations and structure formation can occur [4]. The contribution to entropy production density due to the diffusion of an ideal solution is the product of the diffusive flux and the corresponding driving force:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the diffusivity and concentration of the k-th chemical species; the gas constant R, to be multiplied on the right-hand side, is chosen as 1 for simplicity of calculation [3,7]. Unlike previous studies [18,19], the contribution of entropy flow into and out of the system is not considered.

Irreversible chemical reaction systems are commonly used to model pattern formation, with the reverse reactions often omitted, as mentioned earlier. Accordingly, entropy production diverges for irreversible chemical reactions, as seen in Equation (4). This divergence prevents the calculation and use of entropy production in irreversible reactions [19]. Therefore, reversible chemical reactions, such as the reversible Gray–Scott model, are considered, to overcome this issue. It was first introduced by Mahara et al. in [3,17]. However, Escher, Ross [13], and Yoshida [14] had already considered a similar model for CSTR. The model consists of a cubic autocatalytic reaction $[eqn]$ coupled with a first-order reaction $[eqn]$ corresponding to the following system of equations:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where U, V, and P are the 2D concentrations of the involved chemical species, $[eqn]$ and $[eqn]$ are the rate constants of the forward and backward reactions, respectively, with $[eqn]$ , and f is the feeding rate constant. When backward reaction rates vanish, the original Gray–Scott model is recovered. The second reaction rate and the feeding rate describe the system by fixing the first reaction rate, the diffusion coefficients, and the backward rates. This system can produce many different patterns, such as spots and lines, and there are even indications of regions where different patterns can coexist [15]. The stability analysis of the system and the regions of the phase diagram $[eqn]$ for which this system can form patterns in one and two dimensions are well-documented [15,17,20]. Studies of the phase diagram of the system indicate that it is monostable outside the saddle-node bifurcation curve but bi-stable between the saddle-node bifurcation curve and the Hopf bifurcation curve. Heterogeneous structures are often found precisely around these bifurcation curves [20].

For simplicity and to minimize the system parameters when solving the equations, the control parameters chosen were $[eqn]$ for the first direct reaction, $[eqn]$ for the reverse reactions, and the diffusion coefficients $[eqn]$ , $[eqn]$ and $[eqn]$ . The initial condition was taken as $[eqn]$ , as adopted by previous contributions [3,20]. The system size in 2D is $[eqn]$ with 2500 time steps, and Neumann boundary conditions were used.

The entropy production density for this 2D system is the sum of these two terms:

$$[eqn]$$ $$[eqn]$$

whereas total entropy production is

$$[eqn]$$

In reference [17], it was observed that entropy production from reactive terms is approximately 20 times larger than that arising from the diffusion contribution, yet the latter remains important in affecting the dynamics of pattern formation. For this reason, we will only study the diffusive contribution, since the geometry of the boundaries will primarily influence this term. Figure 1 shows a typical two-dimensional pattern of this system and the corresponding diffusive contribution to the entropy density.

3. Pattern Formation in Reaction–Diffusion on Narrow Channels

When chemical species diffuse in a confined space, their transport is driven by stochastic dynamics and the available space for the motion, i.e., entropic effects. In the ideal case, it is possible to model the effect of periodic channel wrinkles on the diffusion of a particle as an overdamped Brownian motion in a one-dimensional periodic potential [52]:

$$[eqn]$$

where $[eqn]$ is the friction coefficient, $[eqn]$ is the L-periodic potential, F is the static external force, $[eqn]$ is the Boltzmann constant, T is the temperature, and $[eqn]$ is a Gaussian noise. In the free case, i.e., $[eqn]$ , the Einstein relation $[eqn]$ is recovered. For the case $[eqn]$ , the effective Lifson–Jackson coefficient for external potential is recovered [53]:

$$[eqn]$$

When considering non-interacting point particles diffusing in 2D, the process can be described by the corresponding Fokker–Planck equation associated to Equation (12), where the Neumann boundary conditions model the confinement.

It is well known that the system is well-approximated by 1D Brownian diffusion in the presence of entropic barriers [33]. In narrow channels characterized by a transverse length that is significantly smaller than the length of the channel axis, the diffusion equation can be projected to one effective dimension. In this case, the marginal concentrations of the k-species are

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , are the functions that determine the geometry of the boundaries and allow us to introduce the so-called width function as $[eqn]$ . This projection assumes small amplitude boundary fluctuations $[eqn]$ and fast equilibration of the transverse degrees of freedom or uniform transverse concentrations. In this case, the 2D equation is replaced by the 1D Fick–Jacobs–Zwanzig operator:

$$[eqn]$$

where $[eqn]$ is the position-dependent diffusion coefficient of the species k, which usually depends on variations of the channel’s width function. According to [48], the reaction kinetic function changes correspondingly, as $[eqn]$ . Also in [48], it was found that dispersion relations and the range of unstable modes will be affected by the geometry of the system and, therefore, the pattern formation. The major contribution was to explore the impact of channel geometry on reaction–diffusion systems and, consequently, on instability conditions for systems that are constrained by the effective one-dimensional projection technique.

However, the Fick–Jacobs equation is local and it is always desirable to have an equation defined in some particular domain. In [54], the authors present a theoretical link based on the macroscopic diffusion coefficient that is dependent on the porosity, tortuosity, and constriction parameters, by averaging the spatial-dependent diffusion coefficient obtained in the Fick–Jacobs (FJ) approximation. The authors compared the total fluxes predicted by the FJ reduction scheme and the effective medium theory with the numerical results arising from the two-dimensional diffusion equation, in order to show the use of the effective diffusion coefficient, given as follows:

$$[eqn]$$

where

$$[eqn]$$

In Equation (16), $[eqn]$ is the averaged marginal concentration; $[eqn]$ is the mean-effective diffusion coefficient, in terms of constriction factor $[eqn]$ , tortuosity $[eqn]$ , and porosity $[eqn]$ , where W is the approximate width of the rectangular domain containing the pore. This coefficient is analogous to the Lifson–Jackson averaged coefficient, applicable only to periodic channels [53]. By contrast, the coefficient in Equation (16) can be used for any pore, as long as its length is much greater than its width.

4. Effect of Channel Geometry on Entropy Production

4.1. Fick–Jacobs Approximation for Reaction–Diffusion Equations and Entropy Production

Accordingly, by constraining the dynamics of the reversible Gray–Scott model to a narrow channel domain, the system of equations becomes as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

In this case, $[eqn]$ are the Equation (14) marginal concentrations, noting also that some reaction rates are modified by the channel geometry through the width function.

The chemical contribution to the entropy production density at the channel is

$$[eqn]$$

whereas for the effective diffusive component, the Fick–Jacobs probability current and the corresponding driving force, respectively, are given as follows [51]:

$$[eqn]$$ $$[eqn]$$

and so the contribution to the entropy is $[eqn]$ ; thereby,

$$[eqn]$$ $$[eqn]$$

Note that if the channel does not vary with the longitudinal coordinate we recover the standard 1D case. Instead, if $[eqn]$ , some terms dependent on the variations of the channel geometry arise, i.e., entropic forces modifying the entropy production density emerge.

4.2. Effective Diffusion Coefficient in Terms of Availability, Tortuosity, and Constriction

Using the approach of effective diffusion coefficients based on the channel’s geometric features, the problem is simplified to a one-dimensional system. The averaged concentration is defined as

$$[eqn]$$

From this perspective, the reaction–diffusion system becomes a 1D system with modified diffusion coefficients:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the effective diffusion coefficient for each species, now dependent on geometric factors such as constriction, porosity, and tortuosity, which vary according to the specific shape of each channel.

The density of entropy production from reaction will be the same as in the one-dimensional case:

$$[eqn]$$

whereas the one corresponding to diffusion flux is as follows:

$$[eqn]$$

such that

$$[eqn]$$

It is worth noting that all these equations apply to the averaged concentration. Therefore, to compare with the one-dimensional marginal concentration, it is necessary to use the approximation $[eqn]$ , where $[eqn]$ represents availability and W is the effective width of the rectangle where the process occurs.

The reversible Gray–Scott system was integrated using the same parameters as in the first section, with $[eqn]$ and $[eqn]$ , over a line of 200 points and 5000 steps, with periodic boundary conditions and an initial condition $[eqn]$ perturbed by sinusoidal oscillations. Figure 2 shows patterns that appear in the one-dimensional system after some time and persist for a considerable period.

On the other hand, the system (25)–(27) was solved under the same conditions but in the constrained channel shown in Figure 3, with a width function $[eqn]$ . According to [54], this channel’s geometric parameters are $[eqn]$ and the corresponding effective width is $[eqn]$ . It was found that the pattern’s maximum decreases due to geometric effects, as shown in Figure 4.

Finally, the density of entropy production by diffusion at the channel’s center, as a time function, is presented. Figure 5 displays both cases for the function $[eqn]$ . Again, in the confined case, production is reduced compared to free diffusion. It is noted that there are two moments when this function changes abruptly, corresponding to changes in pattern formation. The monotonous part relates to partial semi-stationary patterns formed in each interval, and the characteristic pattern type is shown in the inset.

5. Discussion

The study of thermodynamics in the context of pattern formation within reaction–diffusion systems offers a fascinating challenge. This is particularly true when studying systems with reversible reactions, where the formation of dissipative structures can be interpreted as the coexistence and competition of two distinct states with high- and low-entropy production. Furthermore, it has been observed that confinement modifies the characteristics of the patterns in irreversible systems, depending on the geometrical parameters of the confining walls. This work examined the reversible Gray–Scott model within a narrow channel. Using the effective Fick–Jacobs model for channel diffusion, we observed that geometric modifications to the pattern are preserved in this reversible system. Additionally, we showed that the entropy production density resulting from diffusive flux is a compelling indicator of the dynamics during structure formation. Notably, while the constrained geometry of the channel influences the value of this function, its qualitative behavior remains consistent.

Our subsequent work may delve into the parameter space to examine the behavior of entropy production across various dissipative structures, similar to what was done in [20]. However, we will also address how the system’s geometry influences this behavior, through local and effective theoretical approaches. It would also be intriguing to explore the impact of geometrical effects on entropy production in non-reversible systems, as in [7,8,10], particularly on systems where the reaction kinetics are influenced by the position and distance between species, as in [49]. A comprehensive approach could involve expanding the parameter space to include the geometrical aspects of confinement, such as slope or curvature within the local Fick–Jacobs approach, as well as tortuosity, availability, or porosity, and constriction in the context of effective theory [46,54]. Such an investigation would require more demanding numerical analyses, which we aim to undertake in forthcoming studies.