Coupled vs decoupled approaches for PDE/ODE
systems modeling intercellular signaling
Thomas Carraro1 , Elfriede Friedmann, Daniel Gerecht
Institute for Applied Mathematics, Heidelberg University, INF 293/294, 69120 Heidelberg,
Germany
Abstract
We consider PDE/ODE systems for the simulation of intercellular signaling in
multicellular environments. Since the intracellular processes are described by
few ODEs per cell, the PDE part dominates the solving effort. Thus, it is not
clear if commonly used splitting methods can outperform a coupled approach.
Based on a sensitivity analysis, we present a systematic comparison between
coupled and decoupled approaches for this class of problems and show numerical
results.
Keywords: Coupled PDE/ODE systems, Sensitivity analysis, Multilevel
preconditioner, Intercellular signaling
1. Introduction
Cellular signaling has been mathematically described by a variety of different
models mostly relying on large systems of ordinary differential equations (ODE)
[1]. These earlier models were extended by partial differential equations (PDE)
5
to accurately consider the spatial concentration distribution, localized effects
and signal ranges [2, 3, 4, 5]. In these cellular signalling models the PDE part
captures the diffusion of the signaling proteins and the ODE part the chemical
processes in some parts of the domain which can be considered in a specific
context as well mixed (e.g. cell surface, nucleus). The coupling between these
1 Corresponding author
E-mail address: [email protected]
Preprint submitted to Journal of Computational Physics
April 17, 2015
10
parts occurs in linear and nonlinear terms and through Robin-type boundary
conditions. The analysis and numerics of such models can not be investigated
by standard methods and depends strongly on the structure of the equations
[6]. In this paper we focus on the numerical treatment of such coupled systems.
There are mainly two strategies for implicit solvers of coupled systems: (1)
15
nonlinear methods among them the nonlinear multigrid method also called
“full approximation scheme” (FAS) [7, 8], (2) linearization based approaches
(Newton-type). These methods can be used in a combined approach, where for
example a Newton-type method is applied as smoother for a FAS or a linear
(or nonlinear) multigrid method is used as a preconditioner for a Newton-type
20
method. Our intention in this work is not particularly the comparison and
discussion of advantages and disadvantages of these strategies that depend on
many aspects like, e.g. the accuracy of the Jacobian approximation [9]. Moreover, since in the considered coupled PDE/ODE system the linearization is not
a critical point we choose a Newton-type method preconditioned by a linear
25
multigrid method and study the effect of splitting the PDE part from the ODE
part in the linearized system of equations. A splitting solution approach is
often used when restrictions on accuracy can be relaxed in order to allow an
easier numerical treatment of complicated problems. Such an approach makes
it possible to reuse existing solvers and is widely used in numerical methods
30
for coupled systems, see [10, 11, 12, 13, 14, 15]. In case of strongly coupled
equations, this strategy can only be implemented at high computational costs
through very small time steps or a higher number of iterations in the splitting
scheme. We consider PDE/ODE systems for the simulation of intercellular signaling in multicellular environments. Since the ODE part does not lead to a
35
large discretization system like the PDE part, it is not clear if a splitting method
can outperform a coupled approach.
In this context, the scope of our work is to present a systematic comparison between coupled and decoupled approaches for this class of problems. The method
is based on a sensitivity analysis to compute the strength of the coupling. Ad-
40
ditionally, we compare a multigrid method in which the coupling is considered
2
only at the coarsest level to a fully coupled approach.
We focus on the solution of local problems. In our test cases we consider
systems with eight cells. Therefore, this solution process can be used for example
as local solver for nonlinear preconditioner of Newton-type methods [16] or
45
domain decomposition methods [17].
Application. We solve a coupled PDE/ODE system describing Interleukin-2
(IL-2) signaling between different types of T cells. The underlying model has
been derived by Busse et al. [5] to study immune regulation. The resulting
system consists of one PDE for the intercellular area and of three ODEs for
50
each biological cell. The coupling occurs in linear and nonlinear terms of the
equations, as well as in the Robin type boundary conditions. This system is
prototypical for all signaling models for intercellular communication by diffusing
messengers. Other coupled PDE/ODE models for cellular signaling describe
the interaction with multiple diffusing proteins and more detailed intracellular
55
processes [4] or allow for receptor gradients on the cell surface [2, 3].
Busse et al. [5] performed two-dimensional (2D) simulations to analyze the
dynamics and pattern formation of intercellular signaling between 170 cells.
The solution of this model in the three-dimensional (3D) case is much more
involved and requires numerical treatment that goes beyond the method used in
60
the cited work. Therefore, we have developed the numerical method, presented
in this paper, to efficiently solve the model in realistic 3D environments with
up to 5000 cells [18]. The numerical results have shown that 3D effects have an
important influence on the cell interaction. Thus, new insights into the study
of the immune response have been gained.
65
Outline. The paper is organized as follows. In section 2 we give an abstract
description of the model. We present the mathematical formulation and the
functional setting. We discretize the coupled PDE/ODE system by the finite
element method (FEM) in section 3. We use a sensitivity approach in section
4 to analyze the coupling of the PDE/ODE system and in section 5 we present
70
different solving approaches for the coupled system. We present the numerical
3
results exemplary for a particular application and discuss numerical aspects in
section 6.
2. Mathematical Models for Intercellular Signaling
A model for intercellular signaling consists of a PDE equation for the inter75
action between the cells in the intercellular area Ω coupled with ODEs for the
intracellular processes. We denote by Nc the number of cells in Ω and indicate
by Γi the boundary of each cell i for i = 1, . . . , Nc . The outer boundary of Ω is
denoted by Γout .
(a) 8 interacting cells with surfaces Γi ,
different colors for different cell types
(b) Intercellular area Ω
(cutting plane through intercellular area)
Figure 1: Visualization of the computational domain
Depending on the type of intercellular signaling, different nonlinear operators
describe the dynamics in the intercellular area (AΩ ), e.g. degradation, the
dynamics on the cell surfaces (AΓi ) of each cell and the intracellular processes
(Bi ). We denote the solution of the PDE part by u and the vector of solutions
4
of the ODE part by v. A general formulation of the model can be written as
∂t u(t, x) − D∆u(t, x) + AΩ (u(t, x)) = 0
for (t, x) ∈ (0, T ] × Ω,
D∂n u(t, x) − AΓi (u(t, x), vi (t)) = 0
for (t, x) ∈ (0, T ] × Γi ,
D∂n u(t, x) = 0
∂t vi (t) + Bi (˜
ui (t), vi (t)) = 0
(1)
for (t, x) ∈ (0, T ] × Γout ,
for t ∈ (0, T ],
with given diffusion coefficient D and initial values u(0, x) = u0 , v(0) = v 0 . We
denote the average of u on the surface of Γi by u
˜i and by vi the associated ODE
values with this cell.
R
u
˜i (t) =
Γi
u(t, s) ds
|Γi |
(2)
Similar systems arise in different applications, e.g. in closed-loop cardiovascular
80
simulations where the PDE part is coupled with the ODE part on the surfaces
for in- and outflow, see Moghadam et al. [10].
Remark 2.1. To study the dynamical process and validate the model we compute the entire trajectory. Nevertheless, the simulations mostly converge to a
stable steady state. Therefore we consider the implications on a coupled solver
85
for a computation of the steady state in section 6.3.1 as well.
3. Discretization
For a variational formulation we introduce the Hilbert space V p , e.g. V p =
H 1 , for the PDE part of the equation and the vector space V o = Rn , where n
denotes the number of ordinary differential equations in the system. We define
90
the product space V := V p × V o .
We consider the implicit Euler method as time stepping scheme, and discretize spatially the computational domain Ω by continuous finite elements.
Remark 3.1. For practical use one applies higher order time marching schemes,
e.g. Crank-Nicolson would be a proper choice for the PDE part of the system.
5
95
The ODE part could be solved by higher order Runge-Kutta or backward differentiation formulas. To simplify the implementation we choose the unconditionally
stable implicit Euler method for both parts. For the goal of this work, the numerical results regarding the strength of the PDE/ODE coupling are even more
meaningful for higher order schemes as well, see Remark 6.2.
Considering a time step k we use the semi-discretized weak formulation of
the equations (1) to compute (un+1 , v n+1 ) ∈ V in each time step for all ϕ ∈ V p :
(un+1 , ϕ)Ω + kD(∇un+1 , ∇ϕ) + k(AΩ (un+1 ), ϕ)Ω
X
+k
(AΓi (un+1 , vin+1 ), ϕ)Γi = (un , ϕ)Ω ,
i≤NC
X
v n+1 + k
(3)
Bi (˜
un+1
, vin+1 ) = v n .
i
i≤NC
100
To discretize spatially the system (3), we define a grid {Tl }0≤l≤L , which
consists of a subdivision of the domain in quadrilaterals cells K. The diameters
of the cells hK define a mesh parameter h by the piece-wise constant function
h|K = hK . The discrete solution component uh is sought in the finite dimensional space Vhp ⊂ V p . We choose Vhp as the space of Q1 -elements, the space
105
of functions obtained by transformations of bilinear polynomials defined on a
reference regular unit cell.
Only the PDE part needs to be discretized by the FEM, but due to the
coupled system the ODE part of the discretized solution vh ∈ V o depends on
the spatial discretization as well. Then we can write the fully discretized version
of system (3) for all ϕh ∈ Vhp as follows:
(un+1
, ϕh )Ω + kD(∇un+1
, ∇ϕh )k + (AΩ (un+1
), ϕh )Ω
h
h
h
X
n+1
n
+k
(AΓi (˜
un+1
h,i , vh,i ), ϕh )Γi = (uh , ϕh )Ω ,
i≤NC
X
vhn+1 + k
(4)
n+1
n
Bi (˜
un+1
h,i , vh,i ) = vh .
i≤NC
A nonlinear coupled system arises as discrete system (4) and needs to be
6
solved for each step of the time marching scheme. For the resulting discrete
nonlinear coupled system we introduce the shorter notation
Ah (u, v) = fh ,
(40 )
Bh (u, v) = gh .
We use the subscript h also for the operator Bh to indicate its dependence on
the spatial mesh discretization through the coupling with the PDE part. We
omit the subscript h for the solution components u and v to simply the notation
110
in the next sections.
4. Sensitivity Analysis of the Coupled System
A splitting solving scheme allows solving both parts of the system with a
solver at hand, possibly tuned to solve that specific part of the problem. The
coupling can be implemented in an external (to the two solvers) framework
115
allowing an easy implementation of the more complex coupled problem.
Let Sh : v 7→ u and Th : u 7→ v denote the solution operator for the decoupled
PDE part and ODE part of the discretized system of equations (40 )respectively.
The first equation, u = Sh (v), is solved for a given value of v, then the second
equation, v = Th (u), is solved with the resulting value of u and the cycle is
iterated until a given tolerance is reached. This process can also be written as
a composition of the two operators:
un+1 = Sh Th (un ) .
(5)
A fixed point iteration to solve the system (40 ) has a slow convergence rate
(typically only linear) and the number of fixed point iterations depends on the
nature of the coupling and the model parameters. Nevertheless, a splitting linear
solver can be considered advantageous as part of a Newton scheme. Thus, in120
stead of solving the update of the solution using the Jacobian of the full system,
one updates iteratively the two decoupled parts. In the following sections we
7
will present these two solution approaches, based on a (quasi) Newton method,
in which we compare two different strategies to solve the Newton update.
In the rest of this section, we present a sensitivity approach to decide whether
125
a fixed point iteration or the full system update should be used. As we show
later, the choice depends on the actual model parameter and the method gives
a quantitative index that can be used for practical implementations. We use
the well known fact that a fixed point iteration can converge only if a specific
condition on the iteration operator is fulfilled.
Considering the formulation (5) we write the Jacobian of the fixed point
operator as
J=
∂Sh ∂Th
∂v ∂u
(6)
According to the Banach fixed point theorem the following criterion has to
be fulfilled for the convergence of the fixed point iteration
kJk < 1,
(7)
in some norm k · k. A more convenient criterion is the substitution of the norm
with the spectral radius of the matrix J:
|λmax (J)| < 1.
130
(8)
If this condition is fulfilled a simple fixed point approach converges, therefore
this criterion has been used, e.g. in [19], to define whether the coupling of
the system (40 ) is strong. As will be shown in section 6.3.1, the considered
PDE/ODE system for intercellular signaling is strongly coupled and thus a full
coupling is more effective than a decoupled method. Decoupled methods can
135
still converge for strong coupled systems if embedded in a Newton’s scheme but
require a large number of fixed point iterations, as we show in our numerical
results.
We proceed now with the calculation of the largest eigenvalue of the Jacobian
8
J. Therefore we differentiate the discretized operators Ah and Bh and obtain
the sensitivity equations
A0h,u (ˆ
u, vˆ)uδv + A0h,v (ˆ
u, vˆ)δv = 0,
∀δv ∈ V o ,
(9)
0
0
Bh,v
(ˆ
u, vˆ)vδu + Bh,u
(ˆ
u, vˆ)δu = 0,
∀δu ∈ Vhp .
(10)
and
where we have used the notation
uδv :=
∂u
(δv),
∂v
vδu =
∂v
(δu)
∂u
for the sensitivities and we have written A0h,u and A0h,v for the derivatives of Ah
0
0
with respect to u and v and analogously Bh,u
and Bh,v
for the derivatives of
140
Bh . In the decoupled system, uδv indicates the variation of the PDE solution
perturbing the solution of the ODE system and equivalently vδu is the variation
of the ODE system for a perturbation of the PDE system.
Since the sensitivities in the linear solver strongly depend on the used time
stepping scheme, we consider only the sensitivities for a computation of the
145
steady state. Then, the equations (9) are stationary PDEs to be solved for
each component of δv, while the ODE part (10) consists of algebraic equations
solved for each δu. Therefore we compute the sensitivity matrices ∂Sh /∂v as
a N o × N p matrix and ∂Th /∂u as a N p × N o matrix, where N o denotes the
number of ODE equations and N p the dimension of the PDE discretization.
150
Remark 4.1. In the PDE/ODE system, which we present in section 6, the
coupling between the two parts appears only at the boundaries Γi and only with
the first two components of v. Thus, the product ∂Sh /∂v ∂Th /∂u decouples into
a block diagonal matrix consisting of 2 × 2 matrices for each biological cell. In
addition, we only need to calculate the sensitivities (10) for the restriction of
155
δu on the boundaries Γi , which are nonetheless algebraic equations, so that the
9
major costs to calculate the sensitivities are given by the PDE part (9).
For nonlinear systems of equations the sensitivity analysis depends on a given
point of linearization (ˆ
u, vˆ). We compute an approximate numerical solution
of the system (40 ) for characteristic values of the parameters and choose the
160
computed solution as point of linearization.
5. Numerical Schemes
In this section we present different approaches for a solver of a strongly
coupled PDE/ODE system. Their numerical comparison is presented in section
6.
165
5.1. Nonlinear solver
Newton-type methods provide a flexible and reliable framework for nonlinear
problems by solving a series of linear equations. We present both a splitting
and a fully coupled solver of the linearized subproblems.
5.1.1. Fully coupled Newton’s method
To apply Newton’s method we linearize the system and solve in each Newton
step the system:


fh − Ah (un , v n )
,

=

0
0
δv n+1
gh − Bh (un , v n )
Bh,u
(un , v n ) Bh,v
(un , v n )

170
A0h,u (un , v n )
A0h,v (un , v n )

δun+1

(11)
to obtain the Newton updates δun+1 and δv n+1 used to calculate the next
iterates un+1 = un + δv n and v n+1 = v n + δv n+1 .
5.1.2. Decoupled inexact Newton’s method
Next, we consider a splitting solving scheme for the linear systems defined
in each Newton-step. In a typical decoupled scheme the two systems are solved
175
iteratively in separate solvers. In each Newton step the coupled system is solved
by this decoupled scheme until the residual of the system fulfills a certain accuracy or a maximum of iterations is reached.
10
A standard algorithm for a decoupled Newton’s method is shown in Algorithm 1 which iterates until the approximated solution fulfills a prescribed
accuracy (T OLnewton ). In this scheme we solve in each time step m with a
Newton-type method the solution for the next time step (un+1 , v n+1 ) by calculating a few iterations of the decoupled subsystems. For each Newton step n
the decoupled system is solved by the following fixed point iteration


A0h,u (un , v n )
0

 

δui+1
fh − Ah (un , v n )

=

0
0
gh − Bh (un , v n ) − Bh,v
(un , v n )δui
Bh,v
(un , v n )
δv i+1
A0h,v (un , v n )
(12)
until the Newton updates (δui+1 , δv i+1 ) fulfill the linear residual of the system
(11) to an accuracy (T OLiter ).
180
A common approach to accelerate such a solution process is a Quasi-Newton
iteration in which the Jacobian matrix is only approximated up to a certain
accuracy. In this way the costs per Newton iteration are reduced, while the
number of Newton iterations increase. A trade-off between accuracy and total
costs can enable a reduction of computing time with respect to a full Newton
185
method. Such a Quasi-Newton scheme is obtained if a low accuracy (T OLiter )
or a small maximum number of fixed point iterations (M AX iter ) is chosen.
Algorithm 1: Decoupled algorithm: inexact Newton scheme
n=0
repeat
i=0
repeat
compute Newton updates (δui+1 , δv i+1 ) by solving (12)
evaluate the residual resiter of the linear system (11)
i=i+1
until resiter < T OLiter or i = M AX iter
update the iterate un and v n by δun+1 and δv n+1
evaluate the residual resnewton of the nonlinear system (40 )
n=n+1
until resnewton < T OLnewton
11
This decoupled method is compared to the fully coupled Newton method for
different parameters in numerical tests of section 6.3.1.
5.2. Multigrid scheme
190
In this section we introduce a multilevel preconditioner which can cope with
the strong coupling between PDE and ODEs. Such a coupling arises in the
solver of the linear subsystems if the fully coupled Newton method is used
instead of a splitting scheme. Coupled problems are commonly preconditioned
by block preconditioning approaches, e.g. by simple block diagonal methods
195
or a preconditioning of the Schur complement [20]. We will not use a block
preconditioning approach because of the small dimension of the ODE part but
instead set up a coupled preconditioner based on the linear multigrid method.
In fact, it is well known [21, 8] that the most efficient preconditioner for the
PDE block is a multilevel preconditioner. In general the number of iterations
of the preconditioned linear solver is independent of the mesh refinement. We
compare this approach to the previously described decoupled solving scheme
in Figure 2. We consider a hierarchy of meshes {Tl }0≤l≤L , where the index
0 denotes the root mesh, i.e. the coarsest mesh from which all other meshes
are derived by refinement. In this section we use the following notation for the
system matrix of (11)

A0,l
h,u
Kl := 
0
Bh,u
A0,l
h,v
0
Bh,v


(13)
where the index l indicates the grid refinement level. Note, that we do not use
0
the notation with superscript l in the blocks of the ODE part. In fact, Bh,v
does
200
0
not depend on mesh level, while Bh,u
does depend on the mesh level through the
coupling term u
˜h on the cell boundary, but we use the following approximation:
To reduce the computational costs, and simplify the implementation, the coupling ODE/PDE block is calculated at each level with the term u
˜h computed
at the finest level. In this way the whole ODE part does not depend on the
205
refinement level l. We have observed that this modification does not influence
12
considerably the performance of the multilevel algorithm.
The multigrid scheme is used as a preconditioner for a Krylov method applied
to the system matrix Kl . We use a generalized minimal residual (GMRES)
210
method because of the asymmetry of the system matrix, but a different Krylov
method as, e.g., the BiCG or BiCGStab would also be appropriate for our
purpose.
We show numerically in section 6.3.1 that the efficiency of the preconditioner
is independent of the mesh size. The work presented here is based for the PDE
215
part on previous work by Janssen and Kanschat [22], while for the coupling
PDE/ODE we present some new results.
5.2.1. Transfer operators
For the transfer operators we use the following notation
Rll−1 : Vl → Vl−1 (restriction),
Pll−1 : Vl−1 → Vl (prolongation).
(14)
The restriction and prolongation operators act only on the PDE part, i.e. the
finite element discretization, while the ODE part is transferred by the identity
in both directions. The restriction and prolongation for the PDE part are implemented as intergrid transfers induced by the natural embedding of hierarchical
meshes [22]. In matrix notation the restriction of the whole residual is given by
application of the operator

Rl−1
 l
0
0
I

,
(15)
where I ∈ Rn×n denotes the identity matrix of the ODE part. The prolongation
of the whole residual is defined analogously.
13
220
5.2.2. Smoother
In case of strong coupled problems a common strategy for the smoothing
process is to consider the full coupling only at the coarsest level and to smooth
the two parts separately (decoupled) on the finer levels.
Since in our case the ODE part is small in comparison with the PDE one,
225
we expect that the marginally more expensive smoothing of the whole coupled
system at all levels would be efficient given the strong coupling. Therefore, we
compare the two strategies of (1) smoothing the whole system or (2) smoothing
only the PDE part. For this comparison every efficient smoother would be
appropriate, we have chosen the incomplete LU factorization (ILU). The two
230
smoother are denoted S1 and S2 :
S1 : Incomplete factorization (ILU) of the whole matrix Kl (13);
S2 : incomplete factorization (ILU) of the PDE block as part of a Block Gauss
Seidel scheme.
6. Numerical Results
235
In this section, we make a comparison between the different numerical schemes
presented in the previous section. The following computations were performed
using the C++ library deal.II, see Bangerth et al. [23], with the UMFPACK
library applied as direct solver on the coarse grid level [24]. Further implementation details can be found in [25].
240
6.1. Model
Exemplary, we focus on a model for signaling of Interleukin-2 (IL-2) between
T cells in the lymph node presented by Busse et al. [5]. Interleukin concentrations in the intercellular area regulate the type and strength of the immune
response. The model consists of a reaction-diffusion equation describing the
245
distribution of IL-2 between the T cells in the intercellular area Ω coupled with
ODEs for the intracellular processes by a Robin boundary condition.
We use the following notation:
14
• u(t, s) : [0, T ] × Ω → R describes the concentration of IL-2 in the intercellular area.
250
• Ri (t), Ci (t) and Ei (t) : [0, T ] → R describe the number of IL-2 receptors
(IL-2R), built IL-2/IL-2R receptor-complexes and internalized complexes
for each of the simulated T cells. The receptors are distributed homogeneously on the cell surfaces.
The mathematical model consists of a PDE
∂t u(t, x) = D∆u(t, x) − kd u(t, x)
for all (t, x) in (0, T ] × Ω,
D∂n u(t, s) = qi (t, s) − kon Ri (t)u(t, s) + kof f Ci (t)
∂n u(t, x) = 0
for all (t, s) in (0, T ] × Γi ,
for all (t, x) in (0, T ] × Γout ,
(16a)
coupled with three ODEs for each T cell
∂t Ri (t) = wi0 + wi1
Ci (t)3
− kon Ri (t)e
ui (t)
K 3 + Ci (t)3
− kiR Ri (t) + kof f Ci (t) + krec Ei (t)
for all cells i = 1, . . . , Nc ,
∂t Ci (t) = kon Ri (t)e
ui (t) − (kof f + kiB )Ci (t),
∂t Ei (t) = kiB Ci (t) − (krec + kdeg )Ei (t),
R
u(t, s) ds
,
u
˜i (t) = Γi
|Γi |
(16b)
for given initial conditions u(0), Ri (0), Ci (0) and Ei (0) for all cells i ≤ Nc . The
255
used parameters are described in Table 4. Details about the biological processes,
parameters and initial values are presented in [5].
We consider two cell types which share the same receptor dynamics but differ
in the IL-2 secretion rate:
• Secreting T cells, which omit IL-2 with the secretion rate qi = 2500 mol/h,
260
• Responding T cells with qi = 0.
15
Remark 6.1. We evaluate this system of equations with the sensitivity analysis
presented in section 4. For different numbers, size and distribution of biological
cells, as well as moderate secretion rates q in the biological range, see [5], there
exists a unique stationary state, thus it is possible to use a stationary solver.
We obtain maximal eigenvalues λ of the sensitivity matrix such that
5 < λ < 100.
Therefore, this analysis shows quantitatively that the interaction between the
PDE and the ODE part of the system of equations is strong.
We consider a test problem for the numerical computations of eight cells equidistantly distributed in a 3D environment. The configuration is displayed in Figure
265
1: the responding T cells are in white and the secreting T cell is highlighted.
For this test problem a unique stationary state exists and we obtain a maximal
eigenvalue of the sensitivity matrix λ = 8.88 which indicates a strong coupling
between the PDE and the ODE part.
6.2. Multigrid preconditioner
270
We compare the different smoothers S1 and S2 in a series of numerical tests
for a stationary solver of system (16a), such that the comparison is not influenced
by a time stepping scheme. We compute the number of GMRES steps over all
newton steps (Σn) and the average reduction rate (r) of the residual in each
GMRES step. We proceed with the Newton scheme until an accuracy of 10−6 is
275
reached. Since the number of Newton steps only depends on the coupling, the
nonlinearity of the equation and the accuracy of the solver, it is independent
of the grid refinement. Each of the Newton steps is solved for an accuracy of
10−11 .
We refine the grid globally until the finest grid with 885673 degrees of free-
280
dom compared to 367 on the coarse grid for the PDE part. Additionally 24
degrees of freedom of the ODE part are coupled to the PDE part. We set in
both smoothers the number of iterations to three. This parameter does not
16
influences our comparison. Next we apply the two smoothers S1 and S2 in the
multilevel scheme. Hence, we compare two main approaches: smoothing only
the PDE part or smoothing the whole matrix. Since the coupling between the
Table 1: Reduction rates of different preconditioners
MG-S1 (ILU)
MG-S2 (Bl.-ILU)
L log10 r
Σs
log10 r
Σs
2
2.00
46
1.41
69
3
1.92
51
1.27
77
4
1.85
54
1.21
81
5
1.81
54
1.15
84
Notation: Σs GMRES iterations
in all Newton steps
r average reduction rate
L refinement level
285
PDE and the ODE part in the system is strong, we expect that a good smoother
acting on the coupled system works better. In fact, it can be observed that the
coupled ILU smoother S1 is 35% more effective than S2 with much higher reduction rates. Therefore we choose to apply the smoother S1 in the following
290
sections.
6.3. Comparison of coupled and decoupled schemes
In this section we compare the described coupled and decoupled scheme in
different test cases. Hence the two approaches are:
• a Newton-type method in which the linearized coupled system (11) is
295
solved by a GMRES solver preconditioned by the multigrid method described in section 5.2 with smoother S1 ,
• a Newton-type method in which the linearized system is solved in a decoupled manner by a fixed point iteration defined by the system (12). The
PDE block is solved by a GMRES solver preconditioned by a multigrid
300
method following the implementation of [22]. The solver of choice for
the symmetric part would be the conjugate gradient method (CG), but
17
we use instead the GMRES method for a direct comparison. Nevertheless, we have observed that, in combination with the preconditioner, both
solvers have similar performance.
305
To make the schemes comparable we use a Newton-type method with the same
accuracy of the GMRES solver of T OLiter = 10−11 for both linearized systems.
In this way the number of Newton steps to solve the nonlinear problem is independent of the approach and we can compare the total number of GMRES
steps to solve for an accuracy of T OLnewton = 10−6 .
Figure 2: Schematic representation of the considered coupled and uncoupled schemes
(a) coupled scheme
(b) decoupled scheme
Newton-method
Newton-method
↓
↓
Krylov-solver
Fixed point iteration
↓
.
Multigrid-preconditioner
Direct solver
&
Krylov-solver
↓
↓
Smoother S1
Multigrid-preconditioner
↓
PDE Smoother
310
6.3.1. Stationary solver
We compare the two solvers to compute the steady state solution of the
system (16a). The tests comprehend simulations with biological parameters
that correspond to strong coupling and simulations with artificial parameters
which correspond to a weakly coupled system.
315
• In the simulations with biological parameters the maximal value of the
eigenvalues of the sensitivity matrix is λ = 8.88. We expect thus a strong
PDE/ODE coupling and hence that the decoupled approach is far less
effective than the coupled one.
18
• A weakly coupled test case is created artificially by increasing the parame320
ter kd from 0.1 to 1000. The consequent increment of the degradation of u
diminishes the influence of the uptake of the cells which depends directly
on the components of v. Thus the PDE part is ‘decoupled’ from the ODE
part: the sensitivity analysis yields a maximal eigenvalue of λ = 0.01,
which indicates that the coupling is very weak.
325
In Table 2 we compare the number of Newton steps (n) and the number of total GMRES iterations (Σs) needed to obtain a solution of accuracy
(T OLnewton ). In each Newton step the decoupled scheme described in Algorithm 1 is iterated until a residual resiter < T OLiter is reached without a
given maximum for the number of fixed point iterations M AX iter . The aver-
330
age GMRES iterations per Newton step is denoted by s¯ and the sum over all
GMRES steps by Σs. We globally refine the coarse grid three times up to a
number of 114929 degrees of freedom.
Table 2: Coupled vs decoupled solver
biological problem λ = 8.88
L
n
decoupled
s¯
Σs
n
coupled
s¯
Σs
modified problem λ = 0.01
n
decoupled
s¯
Σs
coupled
n s¯ Σs
2 7 748 5236 7 6.6 46 3
11
33 3 7
3 7 921 6444 7 7.3 51 3
11
33 3 7
4 7 957 6699 7 7.7 54 3 11.7 35 3 7
Notation: Σs GMRES iterations during all Newton steps
n Newton steps
s¯ average GMRES iterations per Newton step
L refinement level
λ largest eigenvalue of the sensitivity matrices
21
21
21
The results show that the number of Newton steps is independent of the
used solving approach due to the high accuracy of the Jacobian, which in the
335
decoupled method is obtained using a small T OLiter . This accuracy comes at
great cost for the decoupled solver, the sum of GMRES iterations is significantly higher (by a factor of more than 100) for the strong coupled biological
problem. The coupled solving scheme is very efficient both for the strong and
19
weak coupled problem. The multigrid preconditioner of the coupled scheme
340
reduces the number of GMRES iterations even in the strong coupled problem
to around seven. The decoupled approach is more competitive for the weak
coupled problem, though the coupled solver is still 40% faster.
6.3.2. Non-stationary solver
Instead of applying a stationary solver we solve in this section the time-
345
dependent system (4) until a stationary solution is reached (20 hours using a
dimensional time variable). In each time step a coupled non-linear system has
to be solved. In contrast to the stationary case, the strength of the coupling is
reduced when using small time steps.
We choose the maximum number of fixed point iterations in each Newton
350
step (M AX iter ) between one and four and compare the fully coupled Newton
method to the decoupled Quasi-Newton scheme.
Table 3: Decoupled and coupled solving of the non-stationary problem
M AX iter
decoupled
∆t = 0.1
Σn
Σs
∆t = 0.01
Σn
Σs
1
2
3
4
1393 3375 5566 14016
754 3792 3395 15217
547 4363 2880 21153
448 4775 2871 23658
coupled
356 1799 2868 11299
Notation: Σn Newton steps in all time steps
Σs Krylow iterations in all Newton steps
M AX iter maximum of iterations
per Newton step
In Table 3 we report the number of computed newton steps (Σn) and the
number of computed GMRES steps (Σs) over all time steps. The results are
listed for computations on a once refined spatial grid (2189 degrees of freedom)
355
with 200 (∆t = 0.1) or 2000 (∆t = 0.01) time steps.
A higher maximal number of fixed point iterations per Newton step increases
the accuracy of the linear solver and thus reduces the number of Newton steps.
The decoupled solving scheme with a maximum of four fixed point iterations re20
sult in near the same number of Newton steps compared to the coupled solving
360
scheme but with more than twice the number of computed GMRES steps. Nevertheless, it can be observed that the number of total GMRES steps decreases
with a reduced number of allowed fixed point iterations per Newton step. Thus
more than one fixed point iteration per Newton step should be avoided, if the
Quasi-Newton method is still converging (this condition cannot be asserted a
365
priori).
As already remarked, the effectiveness of the decoupled solver depends on
the strength of the coupling and thus on the size of the time step. In fact for
time steps ∆t = 0.1 the coupled solver needs around half of the iterations of
the decoupled solver, while for smaller time steps (∆t = 0.01) the iterations of
370
the coupled solver are reduced by 20% compared to the decoupled solver with
M AX iter = 1.
Remark 6.2. The use of a higher order time scheme, e.g. the Crank-Nicolson
scheme, allows for larger time steps, and hence leads to a stronger coupling
during the time integration. Since the coupled solving method is more effective
375
(even for small time steps) in the implicit Euler scheme, it works even better in
a higher order scheme.
7. Conclusion and Outlook
In this paper, we considered a coupled nonlinear system consisting of a
parabolic partial differential equation and many ordinary differential equations,
380
which emerges e.g. in systems biology by modeling intercellular signaling pathways. We presented the numerical treatment for an application in immunology:
the dynamics of cytokine (Interleukine-2) signaling between different types of T
cells in the lymph node, an intercellular signaling pathway which controls the
immune response.
385
The presented methods are valid for a larger class of signaling pathways
involving more diffusing signaling molecules and more chemical interactions be-
21
tween them or for intracellular signaling pathways like in [6] or for other applications modeled by systems of coupled PDE/ODE.
In this paper, we showed in a systematic way a numerical comparison be390
tween coupled and decoupled approaches for a class of models which special
structure (coupled PDE/ODE system) could be effectively exploited. For numerical tests a subproblem of eight T cells was considered.
Moreover, we
used a sensitivity analysis and its numerical implementation to study the coupling/decoupling strategy. This approach, applied to the model for Interleukin-2
395
signaling, indicated that a coupling strategy is better suited for a biologically
relevant range of model parameters. We implemented a solution method based
on a Newton-type solver with a multigrid preconditioner and showed that the
coupling strategy used for all levels of the multigrid is advantageous. Depending
on the time step length, up to around 50% of the computing time is saved by a
400
reduction of the linear solver iterations.
We remark that discretizations consisting of globally refined space and time
grids have been used for the computations in this publication. As final comment,
we indicate a possible strategy for reducing the computation time additionally
by using local mesh refinement both in space and time. Different time grids for
405
the PDE and the ODE part allow, depending on the strength of the coupling,
to decrease the number of time steps for the computationally expensive PDE
part of the model. The essential question is how to choose the two time grids
without decreasing the overall convergence rate of the method. An a posteriori
error estimator for the errors of the PDE and ODE discretization is necessary
410
to reach a certain accuracy efficiently by iterative adaptive refinement. The use
of a refinement strategy based on such an error estimator allows to control the
two time grids separately and obtain an optimal time discretization for both
parts of the system. The complex realization of such a method goes beyond the
scope of this paper and is subject of our current research.
415
22
Parameters
Table 4: Biological Parameters, see [5]
Symbol
Value
Parameter
qi
D
kd
wi0
wi1
K
kon
kof f
kiR
kiC
krec
kdeg
r
d
0 − 22000 mol./cell/h
36000 µm2 /h
0,1/h
150 mol./cell/h
3000 mol./cell/h
1000 mol./cell
111,6 /nM/h
0,83/h
0,64/h
1,7/h
9/h
5/h
5µm
5µm
IL-2 secretion rate
Diffusion coefficient of IL-2
Extracellular IL-2 degradation
Antigen stimulated IL-2 receptor expression rate
Feedback induced IL-2 receptor expression rate
Half-saturation constant of feedback expression
IL-2 association rate constant to IL-2 receptors
IL-2 dissociation rate constant from IL-2 receptors
Internalization rate constant of IL-2 receptors
Internalization rate constant of receptor complexes
Recycling rate constant of IL-2 receptors
Endosomal degradation constant IL-2 receptors
Cell radius
Cell to cell distance
Acknowledgements
The authors gratefully acknowledge Prof. Thomas H¨ofer, DKFZ & BioQuant
Heidelberg, for his support, for fruitful discussions and provisioning this concrete
420
application in immunology.
T.C. was supported by Deutsche Forschungsgemeinschaft (DFG) through the
project CA 633/2-1.
E.F. and D.G. were supported by the Helmholtz Alliance on Systems Biology
(SB Cancer, Submodule V.7) and D.G. additionally by ViroQuant.
425
References
[1] H. A. Kestler, C. Wawra, B. Kracher, M. K¨
uhl, Network modeling of signal
transduction: establishing the global view, Bioessays 30 (11-12) (2008)
1110–1125.
[2] J. Sherratt, P. Maini, W. J¨ager, W. Muller, A receptor based model for
430
pattern formation in hydra, Forma 10 (2) (1995) 77–95.
23
[3] A. Marciniak-Czochra, Receptor-based models with diffusion-driven instability for pattern formation in hydra, Journal of Biological Systems 11 (03)
(2003) 293–324.
[4] E. Friedmann, A. C. Pfeifer, R. Neumann, U. Klingm¨
uller, R. Rannacher,
435
Interaction between experiment, modeling and simulation of spatial aspects
in the JAK2/STAT5 Signaling pathway, Springer, 2013.
[5] D. Busse, M. de la Rosa, K. Hobiger, K. Thurley, M. Flossdorf, A. Scheffold,
T. H¨
ofer, Competing feedback loops shape IL-2 signaling between helper
and regulatory T lymphocytes in cellular microenvironments, Proceedings
440
of the National Academy of Sciences.
[6] E. Friedmann, R. Neumann, R. Rannacher, Well-posedness of a linear
spatio-temporal model of the jak2/stat5 signaling pathway, CMA 15 (2)
(2013) 76–102.
[7] A. Brandt, Multi-level adaptive solutions to boundary-value problems,
445
Mathematics of Computation 31 (138) (1977) 333–390.
[8] W. Hackbusch, Multi-grid methods and applications., Springer Series in
Computational Mathematics, 4. Berlin etc.: Springer- Verlag. XIV, 377,
1985.
[9] D. J. Mavriplis, An assessment of linear versus nonlinear multigrid methods
450
for unstructured mesh solvers, Journal of Computational Physics 175 (1)
(2002) 302–325.
[10] M. E. Moghadam, I. E. Vignon-Clementel, R. Figliola, A. L. Marsden, A
modular numerical method for implicit 0D/3D coupling in cardiovascular
finite element simulations, Journal of Computational Physics 244 (0) (2013)
455
63 – 79, multi-scale Modeling and Simulation of Biological Systems.
[11] C. Farhat, M. Lesoinne, Fast staggered algorithms for the solution of threedimensional nonlinear aeroelastic problems, Numerical Unsteady Aerodynamics and Aeroelastic Simulation (1998) 7–1.
24
[12] C. A. Felippa, K. Park, C. Farhat, Partitioned analysis of coupled me460
chanical systems, Computer methods in applied mechanics and engineering
190 (24) (2001) 3247–3270.
[13] D. Mok, W. Wall, Partitioned analysis schemes for the transient interaction of incompressible flows and nonlinear flexible structures, Trends in
computational structural mechanics (2001) 689–698.
465
[14] M. Heil, Stokes flow in an elastic tubea large-displacement fluid-structure
interaction problem, International journal for numerical methods in fluids
28 (2) (1998) 243–265.
[15] G. Seemann, F. Sachse, M. Karl, D. Weiss, V. Heuveline, O. D¨ossel, Framework for modular, flexible and efficient solving the cardiac bidomain equa-
470
tions using petsc, in: A. D. Fitt, et al. (Eds.), Progress in Industrial Mathematics at ECMI 2008, Mathematics in Industry, Springer Berlin Heidelberg,
2010, pp. 363–369.
[16] X.-C. Cai, D. E. Keyes, L. Marcinkowski, Nonlinear additive Schwarz preconditioners and application in computational fluid dynamics., Int. J. Nu-
475
mer. Methods Fluids 40 (12) (2002) 1463–1470.
[17] A. Quarteroni, A. Valli, Domain decomposition methods for partial differential equations, Tech. rep., Oxford University Press (1999).
[18] K. Thurley, D. Gerecht, E. Friedmann, T. H¨ofer, Three-dimensional gradients of cytokine signaling between t cells, PLoS Computational Biology,
480
accepted.
[19] R. Haftka, J. Sobieszczanski-Sobieski, S. Padula, On options for interdisciplinary analysis and design optimization, Structural optimization 4 (1992)
65–74.
[20] J. Mandel, On block diagonal and schur complement preconditioning, Nu-
485
merische Mathematik 58 (1) (1990) 79–93.
25
[21] J. H. Bramble, Multigrid Methods, Longman Scientific and Technical, London, 1993.
[22] B. Janssen, G. Kanschat, Adaptive multilevel methods with local smoothing for h1 - and hcurl -conforming high order finite element methods, SIAM
490
J. Sci. Comput. 33 (2011) 2095–2114.
[23] W. Bangerth, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, M. Maier,
B. Turcksin, T. D. Young, The deal.II library, version 8.2, Archive of
Numerical Software 3.
[24] T. A. Davis, Algorithm 832: Umfpack v4. 3—an unsymmetric-pattern mul-
495
tifrontal method, ACM Transactions on Mathematical Software (TOMS)
30 (2) (2004) 196–199.
[25] D. Gerecht, Adaptive finite element simulations of coupled PDE/ODE systems modeling intercellular signaling, Ph.D. thesis, University Heidelberg
(2015).
26