2.9 Reaction-Diffusion Model on Large Spatial Scales
2.9.1 Existence of Nontrivial Spatio-Temporal Structures of SolutionsSolutions
To motivate the expectation of the existence of traveling wave solutions, we first consider a reaction-diffusion equation in one spatial dimension
∂tu=∂x2u+f(u), x∈R, t >0, (2.77)
with
f(u) = u(1−u)(u−α),0< α < 1
2. (2.78)
Equation (2.77) together with (2.78) is a simple form of the Fitzhugh-Nagumo equation.
The traveling wave ansatzu(x, t) = φ(x+ct) =:φ(s) leads to the ODE
φ00(s)−cφ0(s) +f(φ(s)) = 0, −∞< s <∞. (2.79) In [29, Section 5.4], the existence of a solution φ(s) to (2.79) for c > 0 with φ(s) → 1 for s → ∞ and φ(s) → 0 for s → −∞ has been proven. This solution corresponds to a heteroclinic orbit between the two steady statesu0 = 0 andu1 = 1, and is given by
φ(s) = 1
1 + exp(−s/√
2) and c=√ 2
1 2 −α
. (2.80)
The two steady statesu0 andu1 are separated by the steady state u=α,0< α <1/2.
Now, the reaction-diffusion system (2.14) established in Section 2.3 possesses a similar structure as (2.77), (2.78) and the existence of a traveling wave solution connecting the two asymptotically stable steady states can be expected.
In order to analyze the spatio-temporal structure of the solution, we solve the reaction-diffusion system (2.75) with homogeneous Neumann boundary conditions (2.76) and the initial conditions
X(r,0) =X(l) for 0≤r <1000 and X(r,0) =X(d) for 1000≤r ≤2000
in the domain [0,2000] numerically. For solving the system, we again use the implicit Euler method together with the Finite Difference method for the spatial, equidistant discretization introduced in Appendix A.4. We set the distance between two nodes of the spatial discretization to ∆r = 0.1 and the size of the time step is ∆t = 50·∆r2 = 0.5.
This choice is possible since the CFL condition does not need to be satisfied for the implicit Euler method. We consider the one-dimensional case and the parameter setting given in (2.8). The numerical solution of the reaction-diffusion system presents a traveling wave solution connecting the two asymptotically stable steady states, see Figure 2.14. We observe that the traveling wave solution moves from the right to the left, i.e., the death state invades the region initially occupied by the life state. Actually, the existence of a traveling wave solution does not depend on the exponent n. For this purpose, we solve the reaction-diffusion system (2.75) for n= 2 and n = 3 numerically, see Figure 2.15 and Figure 2.16. Obviously, the shape of the traveling wave solution depends on the exponent n since the death state concentrations depend on n but the life state concentrations do not depend on n, cf. Section 2.4. Besides the shape, the velocity of the traveling wave solution varies with the exponent n. To measure the velocity, we consider the smallest value for r where the concentration Xa of active caspase 8 takes a value larger than Xamax/2. Every 60 time steps, the coordinate of this point is stored and the comparison of two consecutive coordinates yields the velocity of the traveling wave for each time interval. Since we determine the velocity for different times, we first ascertain that the velocity remains constant except for numerical inaccuracies. Averaging all the computed velocities, the numerical experiments finally give an approximate velocity|¯c| ≈0.2118 for n= 2, |¯c| ≈0.1452 for n = 2.5 and |¯c| ≈ 0.0935 for n = 3. Then, taking L= 10µm, the corresponding physical velocities are|cphys| ≈0.3813µm/s forn= 2,|cphys| ≈0.2614µm/s for n= 2.5 and |cphys| ≈ 0.1684µm/s for n = 3. We see that the velocity decreases with an increasing exponent n and that the physical velocities are of a realistic biological order. In order to study the robustness concerning variations in the parameter values we hold the parameter values ˜kc1 = ˜kc2 = 0.444,˜kp1 = ˜kp2 = 5.556·10−4 and choose four
0 500 1000 1500 2000 0
0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
Figure 2.14: Traveling wave solution which moves from the right to the left for n = 2.5 in cased= 1 for the radially symmetric reaction-diffusion system (2.75) after 2400 (left top), 3600 (right top), 4800 (left bottom) and 6000 (right bottom) time steps. Red line: caspase 8, blue line: caspase 3, green line: pro-caspase 8, black line: pro-caspase 3
arbitrary degradation rates ˜kd1 = 0.0028,˜kd2 = 0.0033,k˜d3 = 0.0022 and ˜kd4 = 0.0028.
Again, the numerical simulation shows ford= 1 the existence of traveling wave solutions connecting the two asymptotically stable steady states. But for this special choice of degradation rates, we observe a non-monotone traveling wave solution, see Figure 2.17.
For n = 2 and n = 3, the traveling wave solution has a similar structure but, due to the difference in the death state concentrations, the traveling wave solutions differ in their shape. Again, we determine the velocity of the traveling waves. We obtain the approximate velocities |¯c| ≈ 0.2347 for n = 2, |¯c| ≈ 0.1658 for n = 2.5 and |¯c| ≈ 0.1113 for n = 3. The corresponding physical velocities are |cphys| ≈ 0.4225µm/s for n = 2,
|cphys| ≈ 0.2985µm/s for n = 2.5 and |cphys| ≈ 0.2004µm/s for n = 3. Thus, for fixed exponentn, the traveling wave solution is faster in case of different degradation rates than in case of equal degradation rates. This is not surprising since the steady states depend on the degradation rates and on the other hand, the steady states influence the velocity of the traveling wave, cf. (2.80).
Remark 2.15. In Section 2.8, we stated Theorem 2.14 with the proposition that the
0 500 1000 1500 2000 0
0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
Figure 2.15: Traveling wave solution which moves from the right to the left for n = 2 in cased= 1 for the radially symmetric reaction-diffusion system (2.75) after 2400 (left) and 6000 (right) time steps. Red line: caspase 8, blue line: caspase 3, green line: pro-caspase 8, black line: pro-caspase 3.
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
Figure 2.16: Traveling wave solution which moves from the right to the left for n = 3 in cased= 1 for the radially symmetric reaction-diffusion system (2.75) after 2400 (left) and 6000 (right) time steps. Red line: caspase 8, blue line: caspase 3, green line: pro-caspase 8, black line: pro-caspase 3.
gradient of the caspase concentrations is balanced by the diffusion on a small time scale.
But now, the eigenvalues of the one-dimensional Laplace operator with homogeneous Neumann boundary conditions on the domain Ω = [0,2000] are λk = k2·π2/20002, especially λ1 = π2/20002 ≈ 2.4674·10−6 and with M = 0.0917 and ˜d = 1, we obtain σ = ˜d·λ1 − M ≈ −0.0917 < 0 and thus, the assumptions for Theorem 2.14 are not fulfilled. So, the numerical result that shows the existence of a traveling wave does not contradict the theoretical results from Section 2.8.
Remark 2.16. In case d= 2 and d = 3, we obtain solutions with a similar structure as the traveling wave solution. We call them front-like solutions. But we expect that the velocity of the solution is not constant anymore due to the curvature in the two- and three-dimensional space.
0 500 1000 1500 2000 0
0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
Figure 2.17: Traveling wave solution which moves from the right to the left forn= 2.5 in case d= 1 for the radially symmetric reaction-diffusion system (2.75) with ˜kd1 >k˜d3 and ˜kd2 >k˜d4
after 2400 (left top), 3600 (right top), 4800 (left bottom) and 6000 (right bottom) time steps.
Red line: caspase 8, blue line: caspase 3, green line: pro-caspase 8, black line: pro-caspase 3.
Remark 2.17. The numerical simulation shows in the one-dimensional case the existence of a traveling wave solution connecting the two asymptotically stable steady states. We notice that the structure of the traveling wave solution is by a factor 100 too large to fit into a cell of size 10µm. Scaling the spatial variable by a factor 100, i.e. settingy=x/100, leads to the scaled second derivative uyy = 1002uxx and so, we obtain
ut= ˜Duxx+f(u) = ˜D/1002uyy+f(u) =D∗uyy+f(u)
with D∗ = ˜D/10000. Thus, in the rescaled equation, the diffusion rate is by a factor 10000 smaller than in the original equation. However, the velocity of the traveling wave then will be so small that execution of cell death would take an unphysiologically long time.
At the end, we switch the region initially occupied by the death state and the region initially filled with the life state, respectively, i.e., we choose the initial conditions
X(r,0) =X(d) for 0≤r <1000 and X(r,0) =X(l) for 1000≤r≤2000.
Again, we observe through numerical simulations a traveling wave solution. But this time, the solution moves from the left to the right, that means, the death state again invades the region initially occupied by the life state, see Figure 2.18.
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
0 500 1000 1500 2000
0 0.05 0.1 0.15 0.2 0.25
r (x 10 um)
caspase concentration (x 100 nM)
Figure 2.18: Traveling wave solution which moves from the left to the right forn= 2.5 in case d= 1 for the radially symmetric reaction-diffusion system (2.75) with ˜kd1 = ˜kd2 = ˜kd3 = ˜kd4 after 2400 (left top), 3600 (right top), 4800 (left bottom) and 6000 (right bottom) time steps.
Red line: caspase 8, blue line: caspase 3, green line: pro-caspase 8, black line: pro-caspase 3.
However, this behavior depends on the special choice of initial condition. For initial conditions of the form
X(r,0) =X(d) for 0≤r < r∗ and X(r,0) =X(l) for r∗ ≤r≤2000,
0< r∗ <2000, the formation of a traveling wave solution is observable if r∗ is larger than a critical value rc. In the next section, we investigate the situation when r∗ < rc and determine the critical value rc, so it exists.