Journal of Nonlinear Science (2021) 31:59
https://doi.org/10.1007/s00332-021-09714-4
Diffusivity Estimation for Activator–Inhibitor Models:
Theory and Application to Intracellular Dynamics of the
Actin Cytoskeleton
Gregor Pasemann1·Sven Flemming2·Sergio Alonso3·
Carsten Beta2·Wilhelm Stannat1
Received: 7 October 2020 / Accepted: 12 April 2021 / Published online: 3 May 2021
© The Author(s) 2021
Abstract
A theory for diffusivity estimation for spatially extended activator–inhibitor dynamics
modeling the evolution of intracellular signaling networks is developed in the math-
ematical framework of stochastic reaction–diffusion systems. In order to account for
model uncertainties, we extend the results for parameter estimation for semilinear
stochastic partial differential equations, as developed in Pasemann and Stannat (Elec-
tron J Stat 14(1):547–579, 2020), to the problem of joint estimation of diffusivity and
parametrized reaction terms. Our theoretical findings are applied to the estimation of
effective diffusivity of signaling components contributing to intracellular dynamics of
the actin cytoskeleton in the model organism Dictyostelium discoideum.
Keywords Parametric drift estimation ·Stochastic reaction–diffusion systems ·
Maximum likelihood estimation ·Actin cytoskeleton dynamics
Communicated by Philipp M Altrock.
BGregor Pasemann
Sven Flemming
Sergio Alonso
Carsten Beta
Wilhelm Stannat
1Institute of Mathematics, Technische Universität Berlin, 10623 Berlin, Germany
2Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam, Germany
3Department of Physics, Universitat Politècnica de Catalunya, 08034 Barcelona, Spain
123
59 Page 2 of 34 Journal of Nonlinear Science (2021) 31 :59
Mathematics Subject Classification 60H15 ·62F10
1 Introduction
The purpose of this paper is to develop the mathematical theory for statistical infer-
ence methods for the parameter estimation of stochastic reaction–diffusion systems
modeling spatially extended signaling networks in cellular systems. Such signaling
networks are one of the central topics in cell biology and biophysics as they provide
the basis for essential processes including cell division, cell differentiation, and cell
motility (Peter 2017). Nonlinearities in these network may cause rich spatiotemporal
behavior including the emergence of oscillations and waves (Beta and Kruse 2017).
Furthermore, alterations and deficiencies in the network topology can explain many
pathologies and play a key role in diseases such as cancer (Condeelis et al. 2005). Here
we present a method to estimate both diffusivity and reaction terms of a stochastic
reaction–diffusion system, given space–time structured data of local concentrations
of signaling components. We mainly focus on the estimation of diffusivity, whose
precision can be increased by simultaneous calibration of the reaction terms.
To test this approach, we use fluorescence microscopy recordings of the actin
dynamics in the cortex of cells of the social amoeba Dictyostelium discoideum,a
well-established model organism for the study of a wide range of actin-dependent
processes (Annesley and Fisher 2009). A recently introduced stochastic reaction–
diffusion model could reproduce many features of the dynamical patterns observed in
the cortex of these cells including excitable and bistable states (Alonso et al. 2018;
Flemming et al. 2020; Moreno et al. 2020). In combination with the experimental data,
this model will serve as a specific test case to exemplify our mathematical approach.
Since in real-world applications the available data will not allow for calibrating and
validating detailed mathematical models, in this paper we will be primarily interested
in minimal models that are still capable of generating all observed dynamical fea-
tures at correct physical magnitudes. The developed estimation techniques should in
practice be as robust as possible w.r.t. uncertainty and even misspecification of the
unknown real dynamics.
The impact of diffusion and reaction in a given model will be of fundamentally
different structure and it is one of the main mathematical challenges to separate
these impacts in the data in order to come to valid parameter estimations. On the
more mathematical side, diffusion corresponds to a second-order partial differential
operator—resulting in a strong spatial coupling in the given data, whereas the reaction
corresponds to a lower order, in fact 0 order, in general resulting in highly nonlinear
local interactions in the data. For introductory purposes, let us assume that our data
is given in terms of a space- and time-continuous field X(t,x)on [0,T]×D, where
Tis the terminal time of our observations and D⊂R2a rectangular domain that
corresponds to a chosen data segment in a given experiment. Although in practice the
given data will be discrete w.r.t. both space and time, we will be interested in applica-
tions where the resolution is high enough in order to approximate the data by such a
continuous field. Our standing assumption is that X(t,x)is generated by a dynamical
123
Journal of Nonlinear Science (2021) 31 :59 Page 3 of 34 59
system of the form
∂tX(t,x)=θ0X(t,x)+FX(t,x), (1)
where is the Laplacian, given by X(t,x)=∂2
x1X(t,x)+∂2
x2X(t,x),x=(x1,x2),
which captures the diffusive spreading in the dynamics of X(t,x). The intensity of
the diffusion is given by the diffusivity θ0. Finally, Fis a generic term, depending on
the solution field X(t,x), which describes all non-diffusive effects present in X(t,x),
whether they are known or unknown. A natural approach to extract θ0from the data
is to use a “cutting-out estimator” of the form
ˆ
θ0=T
0DY(t,x)∂tX(t,x)dxdt
T
0DY(t,x)X(t,x)dxdt
=T
0Y,∂
tXdt
T
0Y,Xdt
,(2)
where Y(t,x)is a suitable test function. In the second fraction of (2), we use the
functional form for readability. In particular, we write X=X(t,x)for the solution
field. We will also write Xt=X(t,·)for the (spatially varying) solution field at a
fixed time t. In order to ease notation, we will use this functional form from now on
throughout the paper. It is possible to derive (2) from a least squares approach by
minimizing θ→∂tX−θX2with a suitably chosen norm. If the non-diffusive
effects described by Fare negligible, we see by plugging (1)into(2) that ˆ
θ0is close
to θ0. If a sound approximation Fto Fis known, the estimator can be made more
precise by substituting ∂tXby ∂tX−FXin (2). A usual choice for Yis a reweighted
spectral cutoff of X, which leads to the spectral approach described below.
Under additional model assumptions, e.g., if (1) is in fact a stochastic partial differ-
ential equation (SPDE) driven by Gaussian white noise, a rather developed parameter
estimation theory for θ0has been established in Pasemann and Stannat (2020)onthe
basis of maximum likelihood estimation (MLE).
In this paper, we are interested in further taking into account also those parts of FX
corresponding to local nonlinear reactions. As a particular example, we will focus on
a recently introduced stochastic reaction–diffusion system of FitzHugh–Nagumo type
that captures many aspects of the dynamical wave patterns observed in the cortex of
motile amoeboid cells (Flemming et al. 2020),
∂tU=DUU+k1U(u0−U)(U−u0a)−k2V+ξ, (3)
∂tV=DVV+(bU −V). (4)
Here, we identify θ0=DUand the only observed data is the activator variable, i.e.,
X=U.
Therefore, in this example the non-diffusive part of the dynamics will be further
decomposed as F=F+ξ, where ξis Gaussian white noise and F=F(U)encodes
the non-Markovian reaction dynamics of the activator. The inhibitor component V
in the above reaction–diffusion system is then incorporated for minimal modeling
purposes to allow the formation of traveling waves in the activator variable Uthat are
123
59 Page 4 of 34 Journal of Nonlinear Science (2021) 31 :59
indeed observed in the time evolution of the actin concentration. This model and its
dynamical features is explained in detail in Sect. 2.1.1.
As noted before, it is desirable to include this additional knowledge into the
estimation procedure (2) by subtracting a suitable approximation Fof the—in
practice—unknown F. Although (3), (4) suggest an explicit parametric form for F,
it is a priori not clear how to quantify the nuisance parameters appearing in the sys-
tem. Thus, an (approximate) model for the data is known qualitatively, based on the
observed dynamics, but not quantitatively. In order to resolve this issue, we extend (2)
and adopt a joint maximum likelihood estimation of θ0and various nuisance parame-
ters.
The field of statistical inference for SPDEs is rapidly growing, see Cialenco (2018)
for a recent survey. The spectral approach to drift estimation was pioneered by Hübner
et al. (1993), Huebner and Rozovskii (1995) and subsequently extended by various
works, see, e.g., Huebner et al. (1997), Lototsky and Rosovskii (1999), Lototsky and
Rozovskii (2000) for the case of non-diagonalizable linear evolution equations. In
Cialenco and Glatt-Holtz (2011), the stochastic Navier–Stokes equations have been
analyzedasafirstexampleofanonlinearevolutionequation.Thishasbeengeneralized
by Pasemann and Stannat (2020) to semilinear SPDEs. Joint parameter estimation for
linear evolution equations is treated in Huebner (1993), Lototsky (2003), see also
Piterbarg and Rozovskii (1996) for a discussion. Besides the spectral approach, other
measurementschemeshavebeenstudied.See,e.g.,PospíšilandTribe(2007),Bibinger
and Trabs (2020), Bibinger and Trabs (2019), Chong (2019a), Chong (2019b), Khalil
and Tudor (2019), Cialenco and Huang (2019), Cialenco et al. (2020), Cialenco and
Kim (2020), Kaino and Uchida (2019) for the case of discrete observations in space
and time. Recently, the local approach has been worked out in Altmeyer and Reiß
(2020) for linear equations, was subsequently generalized in Altmeyer et al. (2020b)
to the semilinear case and applied to a stochastic cell repolarization model in Altmeyer
et al. (2020a).
The paper is structured as follows: In Sect. 2, we give a theory for joint diffusivity
and reaction parameter estimation for a class of semilinear SPDEs and study the spatial
high-frequency asymptotics. Special emphasis is put on the FitzHugh–Nagumo sys-
tem. In Sect. 3, the biophysical context for these models is discussed. The performance
of our method on simulated and real data is evaluated in Sect. 4.
2 Maximum Likelihood Estimation for Activator–Inhibitor Models
In this section, we develop a theory for parameter estimation for a class of semilinear
SPDE using a maximum likelihood ansatz. The application we are aiming at is an
activator–inhibitor model as in Flemming et al. (2020). More precisely, we show
under mild conditions that the diffusivity of such a system can be identified in finite
time given high spatial resolution and observing only the activator component.
123
Journal of Nonlinear Science (2021) 31 :59 Page 5 of 34 59
2.1 The Model and Basic Properties
Let us first introduce the abstract mathematical setting in which we are going to derive
our main theoretical results. We work in spatial dimension d≥1. Given a bounded
domain D=[0,L1]×···×[0,Ld]⊂Rd,L1,...,Ld>0, we consider the following
parameter estimation problem for the semilinear SPDE
dXt=θ0Xt+Fθ1,...,θK(X)dt+BdWt(5)
with periodic boundary conditions for on the Hilbert space H=¯
L2(D)={u∈
L2(D)|Dudx=0}, together with initial condition X0∈H. We allow the non-
linear term Fto depend on additional (nuisance) parameters θ1,...,θKand write
θ=(θ0,...,θK)T,θ1:K=(θ1,...,θK)for short. Without further mentioning it,
we assume that θ∈for a fixed parameter space , e.g., =RK
+.Next,W
is a cylindrical Wiener process modeling Gaussian space–time white noise, that is,
E[˙
W(t,x)]=0 and E[˙
W(t,x)˙
W(s,y)]=δ(t−s)δ(x−y). In order to introduce
spatial correlation, we use a dispersion operator of the form B=σ(−)−γwith
σ>0 and γ>d/4, describing spectral decay of the noise intensity. Here, σis the
overall noise intensity, and γquantifies the decay of the noise for large frequencies in
Fourier space. In addition, γdetermines the spatial smoothness of X, see Sect. 2.1.2.
The condition γ>d/4 ensures that the covariance operator BBTis of trace class,
which is a standard assumption for well-posedness of (5), cf. Liu and Röckner (2015).
Denoteby(λk)k≥0theeigenvaluesof−,orderedincreasingly, with corresponding
eigenfunctions (k)k≥0.Itiswellknown(Weyl1911; Shubin 2001) that λkk2/d
for a constant >0, i.e., limk→∞ λk/(k2/d)=1. The proportionality constant
is known explicitly [see, e.g., Shubin (2001, Proposition 13.1)] and depends on
the domain D.LetPN:H→Hbe the projection onto the span of the first N
eigenfunctions, and set XN:= PNX. For later use, we denote by Ithe identity
operator acting on H.Fors∈R, we write Hs:= D((−)s/2)for the domain of
(−)s/2, which is given by
(−)s/2x=
∞
k=1
λs/2
kk,xk,
and abbreviate |·|
s:= | · |Hsfor the norm on that space whenever convenient. We
assume that the initial condition X0is regular enough, i.e., it satisfies E[|X0|p
s]<∞
for any s≥0, p≥1, without further mentioning it in the forthcoming statements.
We will use the following general class of conditions with s≥0 in order to describe
the regularity of X:
(As)For any p≥1, it holds
Esup
0≤t≤T
|Xt|p
s<∞.(6)
123
Loading more pages...