Building Reaction-Diffusion Models

From Molecular Processes to Spatial Equations

A tissue model needs to represent both movement between locations and reactions within each location. Diffusion redistributes a molecule; production, degradation, and binding change its local amount. A reaction-diffusion model combines these processes in a local balance equation. Its predictions can be uniform across space, or vary with position when transport and local processes interact with nonuniform conditions such as boundary inputs.

The central modeling decision is to define the state variables and specify which processes change each one. The equations below use a freely diffusing ligand and an immobile binding target. They assume a fixed tissue domain, constant diffusion coefficient, and spatially local reactions. Boundary and initial conditions complete the model; their role in defining tissue problems is covered by the earlier material on boundary and initial conditions for tissue models.

For one molecule with constant production rate \(q\) and first-order degradation rate constant \(k_{\mathrm{deg}}\), a basic model is

Here \(c\) has units of amount per volume, \(D\) has units of length squared per time, \(q\) has units of concentration per time, and \(k_{\mathrm{deg}}\) has units of inverse time. Each term therefore has units of concentration per time. Production is positive; first-order degradation is negative and proportional to the concentration present.

Adding Reversible Binding

Suppose the ligand binds reversibly to an immobile target. Let \(b(\mathbf{x},t)\) be bound ligand concentration and \(R_T\) be the local total target concentration, assumed constant. The unbound target concentration is \(R_T-b\). Under mass-action kinetics, association occurs at rate \(k_{\mathrm{on}}c(R_T-b)\), and dissociation occurs at rate \(k_{\mathrm{off}}b\). Both rates are local: molecules bind or dissociate at the position where they are present.

The binding terms have opposite signs in the two equations: association removes free ligand and adds bound ligand, while dissociation does the reverse. The association constant \(k_{\mathrm{on}}\) has units of concentration\(^{-1}\) time\(^{-1}\); \(k_{\mathrm{off}}\) has units of inverse time. This model assumes binding does not itself destroy ligand or target, and that bound ligand does not diffuse. If bound ligand is degraded or mobile, the relevant balance equation must include those processes explicitly.

The target constraint is \(0\leq b\leq R_T\). With nonnegative initial concentrations and physically consistent boundary conditions, the model should preserve nonnegative concentrations and respect that target capacity. These are useful checks on the formulation and on any numerical solution.

Uniform Steady States

A spatially uniform state has no spatial gradients, so \(\nabla^2c=0\). For a no-flux or periodic domain, a uniform initial state remains uniform under spatially uniform parameters. At steady state, both time derivatives are zero. In the model above, setting the bound-state equation to zero gives binding equilibrium:

Define \(K_D=k_{\mathrm{off}}/k_{\mathrm{on}}\), with concentration units. Solving the equilibrium equation for \(b\) gives \(b=R_Tc/(K_D+c)\). At steady state the binding terms in the free-ligand balance also cancel, leaving \(0=q-k_{\mathrm{deg}}c\). For positive \(k_{\mathrm{deg}}\), this gives \(c=q/k_{\mathrm{deg}}\). This result depends on the model assumption that bound ligand is neither degraded nor otherwise removed. Binding changes how ligand is partitioned between free and bound states, but not this steady free concentration under those assumptions.

A Nonuniform Steady Profile

To see how local production and degradation can coexist with spatial variation, consider a one-dimensional tissue layer \(0\leq x\leq L\). Assume no binding, constant production \(q\), first-order degradation \(k_{\mathrm{deg}}c\), and diffusion coefficient \(D\). A maintained boundary concentration at \(x=0\) and a no-flux boundary at \(x=L\) produce the conditions \(c(0)=0\) and \(\partial c/\partial x|_{x=L}=0\). The steady equation is

A spatially uniform steady solution would equal \(q/k_{\mathrm{deg}}\), but it cannot satisfy \(c(0)=0\) when \(q>0\). The boundary condition therefore requires a nonuniform profile. To solve the equation, write \(c=q/k_{\mathrm{deg}}+u\). Substitution gives \(D u''-k_{\mathrm{deg}}u=0\). Its general solution is a linear combination of hyperbolic sine and cosine functions. Applying the two boundary conditions yields

This solution satisfies \(c(0)=0\), because the ratio of the two hyperbolic cosines is then one. Its derivative is zero at \(x=L\), since the derivative of the numerator’s hyperbolic cosine contains \(\sinh(0)=0\). The profile is not uniform: concentration rises away from the cleared boundary and reaches its maximum at the no-flux boundary.

Checks, Alternatives, and Limitations

Before interpreting a reaction-diffusion solution, check that the model and its solution agree in several ways. First, verify the units of every term in each equation. Second, verify reaction signs: a process transferring material between free and bound states must enter those two balances with opposite signs. Third, substitute a proposed steady solution into the differential equation and boundary conditions. Finally, check physical constraints such as nonnegative concentrations and \(b\leq R_T\).

A frequent error is to add a binding sink to the free-ligand equation but omit the corresponding bound-ligand equation. That formulation treats binding as irreversible disappearance even when the process is meant to be reversible. The correction is to include dissociation as a free-ligand source and the matching bound-state balance. Another common error is to assume that a uniform steady state must apply simply because production and degradation are spatially uniform. Uniformity also requires compatible boundary conditions; a fixed zero concentration at one boundary is incompatible with a positive uniform concentration.

A well-mixed compartment model is a useful alternative when spatial gradients are negligible for the question being asked, or when tissue is represented as a single average concentration. It removes the diffusion term and cannot predict position-specific concentrations. A reaction-diffusion model is appropriate when transport and spatially distinct conditions matter, but it requires spatial parameters and boundary conditions. Likewise, explicit binding kinetics are needed when binding transients matter; an equilibrium binding relation is suitable only when binding relaxes sufficiently quickly relative to the changes of interest. Neither choice is universally preferable: the model should retain the mechanisms needed to answer the scientific question without adding unsupported detail.

These equations also simplify real biology. They may assume constant target abundance, homogeneous parameters, linear degradation, and a fixed tissue geometry. Saturable production, target turnover, mobile bound complexes, or cell-specific reactions require additional states or different rate laws. A poor fit or implausible simulation can therefore indicate a numerical issue, a boundary-condition mismatch, or a missing biological process; it is not by itself evidence for which explanation is correct.