Microscopic approach to nonlinear reaction-diffusion: The case of morphogen gradient formation

We develop a microscopic theory for reaction-diffusion (RD) processes based on a generalization of Einstein's master equation [Ann. Phys. 17, 549 (1905)] with a reactive term and show how the mean-field formulation leads to a generalized RD equation with nonclassical solutions. For the $n$th-order annihilation reaction $A+A+A+\cdots +A\to 0$, we obtain a nonlinear reaction-diffusion equation for which we discuss scaling and nonscaling formulations. We find steady states with solutions either exhibiting long-range power-law behavior showing the relative dominance of subdiffusion over reaction effects in constrained systems or, conversely, solutions that go to zero a finite distance from the source, i.e., having finite support of the concentration distribution, describing situations in which diffusion is slow and extinction is fast. Theoretical results are compared with experimental data for morphogen gradient formation.

Figure 1

Steady state: $c\left(r\right)/c\left(0\right)=\text{exp}(-r/r_0)$ [$q=1$, middle (black) curve] and $c\left(r\right)/c\left(0\right)={\{1+[(n-\alpha )/2](r/r_0)\}}^{-2/(n-\alpha )}$ for $n=\alpha +2$ [infinite support, $q=2$, top (red) curve] and $\alpha =n+0.8$ [finite support, $q=0.6$, bottom (blue) curve].

Figure 2

Numerical and analytical solutions for the steady-state profile. (a) Finite support. Numerical solution of the master equation (5), where $p_{j}F=p_j{c}_j^{\alpha -1}$ (with $p_j=0.2$ for $j\in [-2,+2]$ and $\alpha =1.5$) and $G=p_R{c}^n$ (with $n=1$ and $p_R={10}^{-3}$) for $t=50,200,500,800,1.5\times {10}^3,2\times {10}^3,$ and $1.5\times {10}^4$ time steps (symbols), and the analytical steady-state solution (25) [solid (black) curve]. The boundary condition is finite flux at $r=0$ and the initial condition is zero concentration everywhere.(b) Infinite support. Same as (a) except $n=2$ and $p_R={10}^{-2}$ for $t=50,200,500,{10}^3,2\times {10}^3,4\times {10}^3,{10}^4,$and $1\times {10}^5$ time steps (symbols) and the steady-state solution (24) [solid (black) curve]. Note that there are no adjustable parameters in either case.

Figure 3

Numerical and analytical solutions for the finite support profile ($n<\alpha$) with a boundary condition of fixed $c\left(0\right)$. (a) Numerical solution of the master equation (5), where $p_{j}F=p_j{c}_j^{\alpha -1}$ (with $p_j=0.2$ for $j=[-2,+2]$ and $\alpha =1.5)$ and $G=p_R{c}^n$ (with $n=1$ and $p_R={10}^{-3}$) for $t=3\times {10}^2,6\times {10}^2,8\times {10}^2,1.2\times {10}^3,$and $6\times {10}^3$ time steps. (b) Comparison between numerical solution of the master equation (5) for $t=6\times {10}^3$ time steps (open circles) and the analytical steady-state solution (25) (solid curve). Note that there are no adjustable parameters.

Figure 4

Experimental data (black dots) from Han et al. (Fig. 6.A in Ref. [13]) from the fluorescence intensity of the Wg protein (vertical axis, normalized values) versus distance (horizontal axis, in a.u.) measured from the anterior-posterior axis along the dorsoventral direction in the posterior compartment of the Drosophila wild-type wing disc [13]. The solid curve is the best fit of the theoretical steady state (24) with $n-\alpha \simeq 3.8$. For comparison, the dashed curve shows the best-fit exponential profile.

Figure 5

Same as Fig. 4 but for experimental data (black dots) from Han et al. (Fig. 6.B in Ref. [13]) for a mutant strain. The solid curve is the fit of theoretical steady state (24) to the experimental data; because of the obvious asymmetry of the data along the dorsoventral axis, (a) shows a fit based only on the data for negative distances, giving $n-\alpha \simeq 1.7$, and (b) shows a fit to data for positive distances, giving $n-\alpha \simeq 3.3$. In both cases, a best fit to an exponential decay is shown as the dashed curve.