Diffusivity Estimation for Activator-Inhibitor Models: Theory and Application to Intracellular Dynamics of the Actin Cytoskeleton
Gregor Pasemann, Sven Flemming, Sergio Alonso, Carsten Beta, Wilhelm Stannat
DDIFFUSIVITY ESTIMATION FORACTIVATOR-INHIBITOR MODELS: THEORY ANDAPPLICATION TO INTRACELLULAR DYNAMICS OFTHE ACTIN CYTOSKELETON
GREGOR PASEMANN , SVEN FLEMMING , SERGIO ALONSO ,CARSTEN BETA , WILHELM STANNAT Abstract.
A theory for diffusivity estimation for spatially ex-tended activator-inhibitor dynamics modelling the evolution of in-tracellular signaling networks is developed in the mathematicalframework of stochastic reaction-diffusion systems. In order to ac-count for model uncertainties, we extend the results for parameterestimation for semilinear stochastic partial differential equations,as developed in [PS20], to the problem of joint estimation of dif-fusivity and parametrized reaction terms. Our theoretical findingsare applied to the estimation of effective diffusivity of signalingcomponents contributing to intracellular dynamics of the actin cy-toskeleton in the model organism
Dictyostelium discoideum . Introduction
The purpose of this paper is to develop the mathematical theory forstatistical inference methods for the parameter estimation of stochas-tic reaction diffusion systems modelling spatially extended signalingnetworks in cellular systems. Such signaling networks are one of thecentral topics in cell biology and biophysics as they provide the basisfor essential processes including cell division, cell differentiation, andcell motility [DBE + Date : Berlin, May 4, 2020.2010
Mathematics Subject Classification.
Key words and phrases. parametric drift estimation, stochastic reaction-diffusion systems, maximum-likelihood estimation, actin cytoskeleton dynamics. Institute of Mathematics, Technische Universit¨at Berlin, 10623 Berlin, Germany Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam,Germany Department of Physics, Universitat Politcnica de Catalunya, 08034 Barcelona,Spain a r X i v : . [ q - b i o . Q M ] M a y G. PASEMANN, S. FLEMMING, S. ALONSO, C. BETA, W. STANNAT components. We mainly focus on the estimation of diffusivity, whoseprecision can be increased by simultaneous calibration of the reactionterms. To test this approach, we use fluorescence microscopy recordingsof the actin dynamics in the cortex of cells of the social amoeba
Dic-tyostelium discoideum ( D. discoideum ), a well-established model organ-ism for the study of a wide range of actin-dependent processes [AF09].A recently introduced stochastic reaction-diffusion model could repro-duce many features of the dynamical patters observed in the cortex ofthese cells including excitable and bistable states [ASB18, FFAB20].In combination with the experimental data, this model will serve asa specific test case to exemplify our mathematical approach. Since inreal-world applications the available data will not allow for calibrat-ing and validating detailed mathematical models we will be primarilyinterested in this paper in minimal models that are still capable of gen-erating all observed dynamical features at correct physical magnitudes.The developed estimation techniques should be in practice robust asmuch as possible w.r.t. uncertainty and even misspecification of theunknown real dynamics.The impact of diffusion and reaction in a given data set will be of fun-damentally different structure and it is one of the main mathematicalchallenges to separate these impacts in order to come to valid parame-ter estimations. On the more mathematical side, diffusion correspondsto a second order partial differential operator - resulting in a strongspatial coupling in the given data, whereas the reaction corresponds toa lower order, in fact 0-order, in general resulting in highly nonlinearlocal interactions in the data. For introductory purposes, let us as-sume that our data is given in terms of a space and time-continuousfield X ( t, x ) on [0 , T ] × D , where T is the terminal time of our observa-tions and D ⊂ R a rectangular domain that corresponds to a chosendata-segment in a given experiment. Although in practice the givendata will be discrete w.r.t. both space and time, we will be interestedin applications where the resolution is high enough in order to approx-imate the data by such a continuous field. Our standing assumption isthat X is generated by a dynamical system of the form ∂ t X ( t, x ) = θ ∆ X ( t, x ) + F X ( t, x ) , (1)where F is a generic term describing all non-diffusive effects, whetherthey are known or unknown. A natural approach to extract θ fromthe data is to use a ”cutting-out estimator” of the formˆ θ = (cid:82) T (cid:82) D Y ( t, x ) ∂ t X ( t, x )d x d t (cid:82) T (cid:82) D Y ( t, x )∆ X ( t, x )d x d t , (2)where Y is a suitable test function. It is possible to derive this termfrom a least squares ansatz by minimizing θ (cid:55)→ (cid:107) ∂ t X − θ ∆ X (cid:107) with asuitably chosen norm. If the non-diffusive effects described by F are IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 3 negligible, we see by plugging (1) into (2) that ˆ θ is close to θ . Ifa sound approximation F to F is known, the estimator can be mademore precise by substituting ∂ t X by ∂ t X − F X in (2). A usual choicefor Y is a reweighted spectral cutoff of X .Under additional model assumptions, e.g. if (1) is in fact a sto-chastic partial differential equation (SPDE) driven by Gaussian whitenoise, a rather developed parameter estimation theory for θ has beenestablished in [PS20] on the basis of maximum likelihood estimation(MLE).In this paper, we are now interested in further taking into accountalso those parts of F X corresponding to local nonlinear reactions. Asa particular example, we will focus on a recently introduced stochas-tic reaction-diffusion system of FitzHugh-Nagumo type that capturesmany aspects of the dynamical wave patterns observed in the cortex ofmotile amoeboid cells [FFAB20], ∂ t U = D U ∆ U ( t, x ) + k U ( u − U )( U − u a ) − k V + ξ, (3) ∂ t V = D V ∆ V ( t, x ) + (cid:15) ( bU − V ) . (4)Here, we identify θ = D U and the only observed data is the activatorvariable, i.e. X = U .The non-diffusive part of the dynamics will be therefore in this exam-ple further decomposed as F = F + ξ , where ξ is Gaussian white noiseand F = F ( U ) encodes the non-Markovian reaction dynamics of theactivator. The inhibitor component V in the above reaction-diffusionsystem is then incorporated for minimal modelling purposes to allowthe formation of travelling waves in the activator variable U that areindeed observed in the time evolution of the actin concentration.As noted before, it is desirable to include this additional knowledgeinto the estimation procedure (2) by subtracting a suitable approxima-tion F of the - in practice - unknown F . Although (3), (4) suggest anexplicit parametric form for F , it is a priori not clear how to quantifythe nuisance parameters appearing in the system. Thus an (approxi-mate) model for the data is known qualitatively , based on the observeddynamics, but not quantitatively . In order to resolve this issue, we ex-tend (2) and adopt a joint maximum likelihood estimation of θ andvarious nuisance parameters.The field of statistical inference for SPDEs is rapidly growing, see[Cia18] for a recent survey. The spectral approach to drift estima-tion was pioneered by [HKR93, HR95] and subsequently extended byvarious works, see e.g. [HLR97, LR99, LR00] for the case of non-diagonalizable linear evolution equations. In [CGH11], the stochas-tic Navier-Stokes equations have been analyzed as a first example ofa nonlinear evolution equation. This has been generalized by [PS20]to semilinear SPDEs. Joint parameter estimation for linear evolutionequations is treated in [Hue93, Lot03]. Besides the spectral approach, G. PASEMANN, S. FLEMMING, S. ALONSO, C. BETA, W. STANNAT other measurement schemes have been studied. See e.g. [PT07, BT20,BT19, Cho19a, Cho19b, KT19, CH19, CDVK20, CK20, KU19] for thecase of discrete observations in space and time. Recently, the so-calledlocal approach has been worked out in [AR20] for linear equations andwas subsequently generalized in [ACP20] to the semilinear case.The paper is structured as follows: In Section 2 we give a theoryfor joint diffusivity and reaction parameter estimation for a class ofsemilinear SPDEs and study the spatial high-frequency asymptotics.In Section 3, the biophysical context for these models is discussed. Theperformance of our method on simulated and real data is evaluated inSection 4. 2.
Maximum Likelihood Estimation forActivator-Inhibitor Models
In this section we develop a theory for parameter estimation for aclass of semilinear SPDE using a maximum likelihood ansatz. Theapplication we are aiming at is an activator-inhibitor model as in[FFAB20]. More precisely, we show under mild conditions that thediffusivity of such a system can be identified in finite time given highspatial resolution and observing only the activator component.2.1.
The Model and Basic Properties.
Let us first introduce theabstract mathematical setting in which we are going to derive our maintheoretical results. Given a bounded domain D = [0 , L ] ×· · ·× [0 , L d ] ⊂ R d , L , . . . , L d >
0, we consider the following parameter estimationproblem for the semilinear SPDEd X t = ( θ ∆ X t + F θ ,...,θ k ( X )) d t + B d W t (5)with periodic boundary conditions for ∆ on the Hilbert space H =¯ L ( D ) = { u ∈ L ( D ) | (cid:82) D u d x = 0 } , together with initial condition X ∈ H .We write θ = ( θ , . . . , θ K ) T for short and abbreviate θ K = ( θ , . . . , θ K ).Assume that θ ∈ Θ for a suitable parameter space Θ, e.g. Θ = R K + .Here, W is a cylindrical Wiener process, and B = σ ( − ∆) − γ for some γ > d/ σ >
0. Denote by ( λ k ) k ≥ the eigenvalues of − ∆, ordered in-creasingly, with corresponding eigenfunctions (Φ k ) k ≥ . It is well-known[Wey11, Shu01] that λ k (cid:16) Λ k /d for a constant Λ >
0. Let P N : H → H the projection onto the span of the first N eigenfunctions, and set X N := P N X . We write H s := D (( − ∆) s/ ) for the domain of ( − ∆) s/ , s ∈ R , and abbreviate | · | s := | · | H s for the norm on that space when-ever convenient. We assume that the initial condition X is regularenough, i.e. it satisfies E [ | X | ps ] < ∞ for any s ≥ p ≥
1, withoutfurther mentioning it in the forthcoming statements. We will use thefollowing general class of conditions with s ≥ X : IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 5 ( A s ) For any p ≥
1, it holds E (cid:20) sup ≤ t ≤ T | X t | ps (cid:21) < ∞ . (6)Our standing assumption is that there is a unique solution X to (5)such that ( A ) holds. This is a general statement and can be derivede.g. under the assumptions from [LR15, Theorem 3.1.5].2.1.1. An Activator-Inhibitor Model.
An important example for ouranalysis is given by the following system of equations in d ≤ U t = ( D U ∆ U t + k f ( | U t | L , U t ) − k V t )d t + B d W t , (7) d V t = ( D V ∆ V t + (cid:15) ( bU t − V t ))d t, (8)together with sufficiently smooth initial conditions. A detailed dis-cussion of its relevance to cell biology is given in Section 3. Here, f is a bistable third-order polynomial f ( x, u ) = u ( u − u )( u − a ( x ) u ),and a ∈ C b ( R , R ) is a bounded and continuously differentiable functionwith bounded derivative. The boundedness condition for a is not essen-tial to the dynamics of U and can be realized in practice by a suitablecutoff function. For constant a we reobtain the stochastic FitzHugh-Nagumo system.In this particular case, we can find a representation of the abovetype (5) as follows: Using the variation of constants formula, thesolution V to (8) with initial condition V = 0 can be written as V t = (cid:15)b (cid:82) t e ( t − r )( D V ∆ − (cid:15)I ) U r d r . Inserting this representation into (7)yields the following reformulation(9) d U t = (cid:16) D U ∆ U t + k U t ( u − U t )( U t − au ) − k (cid:15)b (cid:90) t e ( t − r )( D V ∆ − (cid:15)I ) U r d r (cid:17) d t + B d W t = ( θ ∆ U t + θ F ( U t ) + θ F ( U t ) + θ F ( U )( t )) d t + B d W t of the activator-inhibitor model (7), (8) by setting θ = D U , θ = k u ¯ a , θ = k , θ = k (cid:15)b , F = 0 for some ¯ a > F ( U ) = − a ( | U | L )¯ a U ( u − U ) , (10) F ( U ) = U ( u − U ) , (11) F ( U )( t ) = − (cid:90) t e ( t − r )( D V ∆ − (cid:15)I ) U r d r. (12) G. PASEMANN, S. FLEMMING, S. ALONSO, C. BETA, W. STANNAT
Here e D V ∆ − (cid:15)I is the semigroup generated by D V ∆ − (cid:15)I . Note that F now depends on the whole trajectory of U , so that the resulting sto-chastic evolution equation (9) is no longer Markovian.For the activator-inhibitor system (7), (8), we can verify well-posednessdirectly. For completeness, we state the optimal regularity results forboth U and V , but our main focus lies on the observed variable X = U . Proposition 1.
Let γ > d/ . Then there is a unique solution ( U, V ) to (7) , (8) . Furthermore, U satisfies ( A s ) for any s < γ − d/ ,and V satisfies ( A s ) for any s < γ − d/ . The proof is deferred to Appendix A.2.1.2.
Basic Regularity Results.
The nonlinear term F is assumed tosatisfy (cf. [ACP20]):( F s,η ) There is b > (cid:15) > | ( − ∆) s − η + (cid:15) F θ K ( Y ) | C (0 ,T ; H ) ≤ c ( θ K )(1 + | ( − ∆) s Y | C (0 ,T ; H ) ) b for Y ∈ C (0 , T ; H s ), where c depends continuously on θ K .In particular, if F ( Y )( t ) = F ( Y t ), this simplifies to | F θ K ( Y ) | s − η + (cid:15) ≤ c ( θ K )(1 + | Y | s ) b (13)for Y ∈ H s . In order to control the regularity of X , we apply a splittingargument (see also [CGH11, PS20, ACP20]) and write X = X + (cid:101) X ,where X is the solution to the linear SPDEd X t = θ ∆ X t d t + B d W t , X = 0 , (14)where W is the same cylindrical Wiener process as in (5), and (cid:101) X solvesa random PDE of the formd (cid:101) X = ( θ ∆ (cid:101) X t + F θ K ( X + (cid:101) X )( t ))d t, X = X . (15) Lemma 2.
The process X is Gaussian, and for any p ≥ , s < γ − d/ : E (cid:20) sup ≤ t ≤ T | X t | ps (cid:21) < ∞ . (16) Proof.
This is classical, see e.g. [DPZ14, LR15]. (cid:3)
Proposition 3. (1)
Let ( A s ) and ( F s,η ) hold. Then for any p ≥ : E (cid:20) sup ≤ t ≤ T | (cid:101) X t | ps + η (cid:21) < ∞ . (17) In particular, if s + η < γ − d/ , then ( A s + η ) is true. IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 7 (2)
Let G : C (0 , T ; H ) ⊃ D ( G ) → C (0 , T ; H ) be any function suchthat ( F s,η ) holds for G . Then for s < γ − d/ and p ≥ : E (cid:20) sup ≤ t ≤ T | G ( X )( t ) | ps + η − (cid:21) < ∞ . (18) In particular, E (cid:90) T | G ( X )( t ) | s + η − d t < ∞ . (19) Proof. (1) For t ∈ [0 , T ] and (cid:15) > | (cid:101) X t | s + η ≤ | S ( t ) X | s + η + (cid:90) t | S ( t − r ) F θ K ( X )( t ) | s + η d r ≤ | X | s + η + (cid:90) t ( t − r ) − (cid:15)/ | F θ K ( X )( t ) | s − η + (cid:15) d r ≤ | X | s + η + 2 (cid:15) T (cid:15) sup ≤ t ≤ T | F θ K ( X )( t ) | s − η + (cid:15) ≤ | X | s + η + 2 (cid:15) T (cid:15) c ( θ K )(1 + | X | C (0 ,T ; H s ) ) b , where θ , . . . , θ K are the true parameters. This implies (17). If s + η < γ − d/ X byLemma 2, and the claim follows.(2) This follows from E (cid:20) sup ≤ t ≤ T | G ( X )( t ) | ps + η − (cid:21) ≤ c E (cid:34)(cid:18) ≤ t ≤ T | X t | s (cid:19) bp (cid:35) < ∞ . (20) (cid:3) Statistical Inference: The General Model.
The projectedprocess P N X induces a measure P θ on C (0 , T ; R N ). Heuristically (see[LS77, Section 7.6.4]), we have the following representation for thedensity with respect to P Nθ for an arbitrary reference parameter θ ∈ Θ:d P Nθ d P Nθ ( X N ) = exp (cid:18) − σ (cid:90) T (cid:10) ( θ − θ )∆ X Nt + P N ( F θ K − F θ K )( X ) , ( − ∆) γ d X Nt (cid:11) + 12 σ (cid:90) T (cid:10) ( θ − θ )∆ X Nt + P N ( F θ K − F θ K )( X ) , ( − ∆) γ (cid:2) ( θ + θ )∆ X Nt + P N ( F θ K + F θ K )( X ) (cid:3)(cid:11) d t (cid:1) . By setting the score (i.e. the gradient with respect to θ of the log-likelihood) to zero, and by formally substituting the (fixed) parameter G. PASEMANN, S. FLEMMING, S. ALONSO, C. BETA, W. STANNAT γ by a (free) parameter α , we get the following maximum likelihoodequations:ˆ θ N (cid:90) T | ( − ∆) α X Nt | H d t − (cid:90) T (cid:104) ( − ∆) α X Nt , P N F ˆ θ N K ( X ) (cid:105) d t = − (cid:90) T (cid:104) ( − ∆) α X Nt , d X Nt (cid:105) , − ˆ θ N (cid:90) T (cid:104) ( − ∆) α X Nt , ∂ θ i P N F ˆ θ N K ( X ) (cid:105) d t + (cid:90) T (cid:104) ( − ∆) α P N F ˆ θ N K ( X ) , ∂ θ i P N F ˆ θ N K ( X ) (cid:105) d t = (cid:90) T (cid:104) ( − ∆) α ∂ θ i P N F ˆ θ N K ( X ) , d X Nt (cid:105) . Any solution (ˆ θ N , . . . , ˆ θ NK ) to these equations is a (joint) maximumlikelihood estimator (MLE) for ( θ , . . . , θ K ). W.l.o.g. we assume thatthe MLE is unique, otherwise fix any solution. We are interested inthe asymptotic behavior of this estimator as N → ∞ , i.e. as more andmore spatial information (for fixed T >
0) is available. While iden-tifiability of θ , . . . , θ K in finite time depends in general on additionalstructural assumptions on F , the diffusivity θ is expected to be identi-fiable in time under mild assumptions. Indeed, the argument is similarto [CGH11, PS20], but we have to take into account the dependenceof ˆ θ N on the other estimators ˆ θ N , . . . , ˆ θ NK . Note that the likelihoodequations give the following useful representation for ˆ θ N :ˆ θ N = − (cid:82) T (cid:104) ( − ∆) α X Nt , d X Nt (cid:105) + (cid:82) T (cid:104) ( − ∆) α X Nt , P N F ˆ θ N K ( X ) (cid:105) d t (cid:82) T | ( − ∆) α X Nt | H d t . (21)Recall that a sequence of estimators (ˆ θ N ) N ∈ N is bounded in probability resp. tight if sup N ∈ N P ( | ˆ θ N | > M ) → M → ∞ . Theorem 4.
Assume that the likelihood equations are solvable for N ≥ N , assume that (ˆ θ Ni ) N ≥ N is bounded in probability for i = 1 , . . . , K .Let α > γ − d/ − / and η, s > such that ( A s ) and ( F s,η ) are truefor any s ≤ s < γ + 1 − d/ . Then the following is true: (1) ˆ θ N is a consistent estimator for θ , i.e. ˆ θ N P −→ θ . (2) If η ≤ d/ , then N r (ˆ θ N − θ ) P −→ for any r < η/d . (3) If η > d/ , then N + d (ˆ θ N − θ ) d −→ N (0 , V ) , (22) with V = 2 θ (4 α − γ + d + 2) / ( T d Λ α − γ +1 (8 α − γ + d + 2)) . IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 9
Lemma 5.
Let α > γ − d/ − / , let further η, s > such that ( A s ) and ( F s,η ) are true for s ≤ s < γ + 1 − d/ . Then (cid:90) T | ( − ∆) α X Nt | H d t (cid:16) E (cid:90) T | ( − ∆) α X Nt | H d t (23) (cid:16) C α N d (2 α − γ +1)+1 (24) in probability, with C α = T Λ α − γ +1 d θ (4 α − γ + 2 + d ) . (25) Proof.
Using Proposition 3 (i), the proof is exactly as in [PS20, Propo-sition4.6]. (cid:3)
Proof of Theorem 4.
Using (21), the proof works exactly as in [PS20].Denote by ˆ θ full ,N the estimator which is given by (21) if the ˆ θ N , . . . , ˆ θ NK are substituted by the true values θ , . . . , θ K . By plugging in the dy-namics (5), we can writeˆ θ full ,N − θ = − c N (cid:82) T (cid:104) ( − ∆) α − γ X Nt , d W Nt (cid:105) (cid:113)(cid:82) T | ( − ∆) α − γ X Nt | H d t . (26)with c N = (cid:113)(cid:82) T | ( − ∆) α − γ X Nt | H d t (cid:82) T | ( − ∆) α X Nt | H d t . By Lemma 5, the rescaled prefactor c N N / /d converges in probabilityto (cid:112) C α − γ /C α . The second factor converges in distribution to a stan-dard normal distribution N (0 ,
1) by the central limit theorem for localmartingales (see [LS89, Theorem 5.5.4 (I)], [JS03, Theorem VIII.4.17]).This proves (22) for ˆ θ full ,N . To conclude, we bound the bias term de-pending on ˆ θ N , . . . , ˆ θ NK as follows, using | P N Y | s ≤ λ ( s − s ) / N | P N Y | s for s < s : Let δ >
0. Then (cid:90) T (cid:104) ( − ∆) α X Nt , P N F ˆ θ N K ( X ) (cid:105) d t ≤ (cid:18)(cid:90) T | ( − ∆) α X Nt | H d t (cid:19) (cid:18)(cid:90) T | ( − ∆) α P N F ˆ θ N K ( X ) | H d t (cid:19) (cid:46) N d (2 α − γ +1)+ (cid:18)(cid:90) T | ( − ∆) α P N F ˆ θ N K ( X ) | H d t (cid:19) (cid:46) N d (2 α − γ +1)+1 − η − δd (cid:18)(cid:90) T | ( − ∆) γ + − d − η − δ P N F ˆ θ N K ( X ) | H d t (cid:19) , so using ( F γ +1 − d/ − δ,η ), (cid:82) T (cid:104) ( − ∆) α X Nt , P N F ˆ θ N K ( X ) (cid:105) d t (cid:82) T | ( − ∆) α X Nt | H d t (cid:46) c (ˆ θ N K ) N − ( η − δ ) /d As c (ˆ θ N K ) is bounded in probability and δ > θ , . . . , θ K is similar. This concludes the proof. (cid:3) It is clear that a Lipschitz condition on F with respect to θ K allowsto bound ˆ θ N − ˆ θ full ,N in terms of | ˆ θ N K − θ K | N − ( η − δ ) /d for δ >
0, usingthe notation from the previous proof. In this case, consistency of ˆ θ Ni , i = 1 , . . . , K , may improve the rate of convergence of ˆ θ N . However,as noted before, in general we cannot expect ˆ θ Ni , i = 1 , . . . , K , to beconsistent as N → ∞ .2.3. Statistical Inference: The Linear Model.
We put particularemphasis on the case that the nonlinearity F depends linearly on itsparameters:d X t = (cid:32) θ ∆ X t + K (cid:88) i =0 θ i F i ( X ) + F ( X ) (cid:33) d t + B d W t . (27)This model includes the FitzHugh-Nagumo system in the form (9).We state an additional verifiable condition, depending on the contrastparameter α ∈ R , which guarantees that the likelihood equations arewell-posed, amongst others.( L α ) The terms F ( Y ) , . . . , F K ( Y ) are well-defined and linearly inde-pendent in L (0 , T ; H α ) for every non-constant Y ∈ C (0 , T ; C ( D )).In particular, condition ( L α ) implies for i = 1 , . . . , K that (cid:90) T | ( − ∆) α F i ( X ) | H d t > . (28)The maximum likelihood equations for the linear model (27) simplifyto A N ( X )ˆ θ N ( X ) = b N ( X ) , (29)where A N ( X ) , = (cid:90) T | ( − ∆) α X Nt | H d t, (30) A N ( X ) ,i = A N ( X ) i, = − (cid:90) T (cid:104) ( − ∆) α X Nt , P N F i ( X ) (cid:105) d t, (31) A N ( X ) i,j = (cid:90) T (cid:104) ( − ∆) α P N F i ( X ) , P N F j ( X ) (cid:105) d t (32) IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 11 for i, j = 1 , . . . , K , and b N ( X ) = − (cid:90) T (cid:104) ( − ∆) α X Nt , d X Nt (cid:105) + (cid:90) T (cid:104) ( − ∆) α X Nt , P N F ( X ) (cid:105) d t, (33) b N ( X ) i = (cid:90) T (cid:104) ( − ∆) α P N F i ( X ) , d X Nt (cid:105) − (cid:90) T (cid:104) ( − ∆) α P N F i ( X ) , P N F ( X ) (cid:105) d t (34)for i = 1 , . . . , K .In order to apply Theorem 4, we have to prove that the estimatorsˆ θ N , . . . , ˆ θ NK are bounded in probability. Thus, the goal of this sectionis to prove: Proposition 6.
In the setting of this section, let η, s > such that ( A s ) and ( F s,η ) are true for s ≤ s < γ + 1 − d/ . For γ − d/ − / <α ≤ ( γ − d/ − / η/ ∧ ( γ − d/ − / η/ , let ( L α ) be true.Then the ˆ θ Ni , i = 0 , . . . , K , are bounded in probability. We note that the upper bound on α can be relaxed in general, de-pending on the exact asymptotic behaviour of A N ( X ) i,i , i = 1 , . . . , K .Proposition 6 together with Theorem 4 gives conditions for ˆ θ N to beconsistent and asymptotically normal in the linear model (27). In par-ticular, we immediately get for the activator-inhibitor model (7), (8),as the linear independency condition ( L α ) is trivially satisfied and η can be chosen arbitrarily close to 2: Theorem 7.
Let γ > d/ . Then ˆ θ N has the following properties in theactivator-inhibitor model (7) , (8) : (1) In d = 1 , let γ − / < α < γ + 1 / . Then ˆ θ N is a consistentestimator for θ , which is asymptotically normal as in (22) . (2) In d = 2 , let γ − < α < γ . Then ˆ θ N is a consistent estimatorfor θ with optimal convergence rate, i.e. N r (ˆ θ N − θ ) P −→ forany r < . The rest of this section is devoted to proving Proposition 6. Thenext statement characterizes the rate of det( A N ( X )). While thereis a trival upper bound det( A N ( X )) (cid:46) a N f N . . . f NK obtained by theCauchy-Schwarz inequality, the corresponding lower bound requiresmore work. Our argument is geometric in nature: A Gramian ma-trix measures the (squared) volume of a parallelepiped spanned by thevectors defining the matrix. The argument takes place in the separableHilbert space L (0 , T ; H α ), and we find a (uniform in N ) lower boundon the ( K + 1)-dimensional volume of the parallelepiped spanned by P N ( − ∆) X, P N F ( X ) , . . . , P N F K ( X ). While the latter K components(and thus their K -dimensional volume) converge, the first component gets eventually orthogonal to the others as N grows. This is formalizedas follows: Lemma 8.
Let η, s > such that ( A s ) and ( F s,η ) are true for s ≤ s < γ + 1 − d/ . Let γ − d/ − / < α ≤ γ − d/ − / η/ . Underassumption ( L α ) , we have det( A N ( X )) (cid:38) a N f N . . . f NK . (35) In particular, A N ( X ) is invertible if N is large enough. A similar argument can be found in [Hue93, Proposition 3.1] forlinear SPDEs.
Proof.
The argument is pathwise, so we fix a realization of X . We ab-breviate (cid:104)· , ·(cid:105) α := (cid:104)· , ·(cid:105) L (0 ,T ; H α ) and (cid:107)·(cid:107) α := | · | L (0 ,T ; H α ) , further wewrite F ( X ) = ∆ X in order to unify notation. By a simple normaliza-tion procedure, it suffices to show thatlim inf N →∞ det (cid:32)(cid:18)(cid:28) P N F i ( X ) (cid:107) P N F i ( X ) (cid:107) α , P N F j ( X ) (cid:107) P N F j ( X ) (cid:107) α (cid:29) α (cid:19) i,j =0 ,...K (cid:33) > . (36)Let (cid:15) > N , N ∈ N with N ≤ N such that: • (cid:107) P N F i ( X ) − F i ( X ) (cid:107) α < (cid:15) (cid:107) F i ( X ) (cid:107) α and (cid:107) F i ( X ) (cid:107) α / (cid:107) P N F i ( X ) (cid:107) α < i = 1 , . . . , K and N ≥ N . This is possible due to α ≤ γ − d/ − / η/ . • (cid:107) P N F ( X ) / (cid:107) P N F ( X ) (cid:107) α (cid:107) α < (cid:15) for N ≥ N . This is possibleas (cid:107) P N F ( X ) (cid:107) α diverges due to α > γ − d/ − / (cid:12)(cid:12)(cid:12)(cid:12)(cid:28) P N F ( X ) (cid:107) P N F ( X ) (cid:107) α , P N F i ( X ) (cid:107) P N F i ( X ) (cid:107) α (cid:29) α (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12)(cid:28) P N F ( X ) (cid:107) P N F ( X ) (cid:107) α , P N F i ( X ) (cid:107) P N F i ( X ) (cid:107) α (cid:29) α (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:28) ( P N − P N ) F ( X ) (cid:107) P N F ( X ) (cid:107) α , ( P N − P N ) F i ( X ) (cid:107) P N F i ( X ) (cid:107) α (cid:29) α (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:13)(cid:13)(cid:13)(cid:13) P N F ( X ) (cid:107) P N F ( X ) (cid:107) α (cid:13)(cid:13)(cid:13)(cid:13) α + (cid:13)(cid:13)(cid:13)(cid:13) ( I − P N ) F i ( X ) (cid:107) P N F i ( X ) (cid:107) α (cid:13)(cid:13)(cid:13)(cid:13) α ≤ (cid:15) + 2 (cid:15) for N ≥ N with i = 1 , . . . , K , where we used the Cauchy-Schwarzinequality and the choice of N , N . Consequently, in order for (36) tobe true, it suffices to show for the (0 , N →∞ det (cid:32)(cid:18)(cid:28) P N F i ( X ) (cid:107) P N F i ( X ) (cid:107) α , P N F j ( X ) (cid:107) P N F j ( X ) (cid:107) α (cid:29)(cid:19) i,j =1 ,...K (cid:33) > . (37) IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 13
By ( L α ), we havedet (cid:32)(cid:18)(cid:28) F i ( X ) (cid:107) F i ( X ) (cid:107) α , F j ( X ) (cid:107) F j ( X ) (cid:107) α (cid:29) α (cid:19) i,j =1 ,...K (cid:33) > , (38)so (37) holds true by continuity for i = 1 , . . . , K . (cid:3) Proof of Proposition 6.
We formally write F ( X ) = ∆ X in order tounify notation. By plugging in the dynamics of X , we see that for i = 0 , . . . , K , b N ( X ) i = (cid:90) T (cid:104) ( − ∆) α − γ P N F i ( X ) , d W Nt (cid:105) + K (cid:88) k =0 θ k A N ( X ) i,k . (39)Thus, with ¯ b N ( X ) i = (cid:82) T (cid:104) ( − ∆) α − γ P N F i ( X ) , d W Nt (cid:105) and θ = ( θ , . . . , θ K ) T ,this read as A N ( X )ˆ θ N ( X ) = ¯ b N ( X ) + A N ( X ) θ, (40)i.e. ˆ θ N − θ = A − N ( X )¯ b N ( X ) . (41)Using the explicit representation for the inverse matrix A N ( X ) − , wehave ˆ θ Nj − θ j = 1det( A N ( X )) K (cid:88) i =0 ( − i + j det( A ( i,j ) N ( X ))¯ b N ( X ) i , (42)where A ( i,j ) N ( X ) results from A N ( X ) by erasing the i th row and the j thcolumn. We abbreviate a Ni := A N ( X ) i,i for i = 0 , . . . , K . Rememberthat by Lemma 5, a N (cid:16) C α N d (2 α − γ +1)+1 . Further, due to α ≤ γ − d/ − / η/ a Ni /a N → i = 1 , . . . , K . Now, bythe Cauchy-Schwarz inequality, each summand in the numerator canbe bounded by terms of the form K (cid:89) k =0 a Nk | ¯ b N ( X ) i | (cid:113) a Ni a Nj (43)for i, j = 0 , . . . , K . Thus, by Lemma 8, in order to prove the claim it re-mains to find a uniform bound in probability for the terms | ¯ b N ( X ) i | / (cid:113) a Ni a Nj .Now, E (cid:2) ¯ b N ( X ) i (cid:3) = E (cid:34)(cid:18)(cid:90) T (cid:104) ( − ∆) α − γ P N F i ( X ) , d W Nt (cid:105) (cid:19) (cid:35) (44) = E (cid:20)(cid:90) T | ( − ∆) α − γ P N F i ( X ) | H d t (cid:21) < ∞ , (45) as α < γ − d/ − / η/
4. Further, ( a N ) − / converges to zero inprobability, while ( a Ni ) − / converges to a positive constant for i =1 , . . . , K . The claim follows. (cid:3) Application to Activator-Inhibitor Models of ActinDynamics
The actin cytoskeleton is a dense polymer meshwork at the inner faceof the plasma membrane that determines the shape and mechanical sta-bility of a cell. Due to the continuous polymerization and depolymer-ization of the actin filaments, it displays a dynamic network structurethat generates complex spatiotemporal patterns. These patterns arethe basis of many essential cellular functions, such as endocytic pro-cesses, cell shape changes, and cell motility [BBPSP14]. The dynamicsof the actin cytoskeleton is controlled and guided by upstream signal-ing pathways, which are known to display typical features of nonequi-librium systems, such as oscillatory instabilities and the emergence oftraveling wave patterns [DBE +
17, BK17]. Here we use giant cells of thesocial amoeba
D. discoideum that allow us to observe these cytoskele-tal patters over larger spatial domains [GEW + + D. dis-coideum in the wave-forming regime for comparison. Images wererecorded by confocal laser scanning microscopy and display the dis-tribution of mRFP-LimE∆, a fluorescent marker for filamentous actin,in the cortex at the substrate-attached bottom membrane of the cell.As individual actin filaments are not resolved by this method, the in-tensity of the fluorescence signal reflects the local cortical density offilamentous actin. Rectangular subsections of the inner part of thecortex of giant cells as displayed in panel (C) were used for data anal-ysis in Section 4.Many aspects of subcellular dynamical patterns have been addressedby reaction-diffusion models. While some models rely on detailed mod-ular approaches [BAB08, DBE + +
10] or Ras signaling [FMU19]. To de-scribe wave patterns in the actin cortex of giant
D. discoideum cells, thenoisy FitzHugh-Nagumo type reaction-diffusion system (3), (4), com-bined with a dynamic phase field, has been recently proposed [FFAB20].
IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 15 (A) Single cell (C) Section of a giant cell(B) Giant cell
ExperimentsSimulations (E) Giant cell (F) Periodic boundaries(D) Single cell
Figure 1.
Actin waves in experiments (top) and modelsimulations (bottom). (A) Normal-sized cell with a cir-cular actin wave. (B) Giant cell with several fragmentedactin waves. (C) Subsection of the cortical area of thegiant cell shown in (B), indicated as dotted green rec-tangle. Experimental images are confocal microscopyrecordings of mRFP-LimE∆ expressing
D. discoideum
AX2 cells, see [GEW + µ m) Details on the numerical implemen-tation can be found in Appendix B.In contrast to the more detailed biochemical models, the structure ofthis model is rather simple. Waves are generated by noisy bistable/exci-table kinetics with an additional control of the total amount of activator U . This constraint dynamically regulates the amount of U around aconstant level in agreement with the corresponding biological restric-tions. Elevated levels of the activator represent typical cell front mark-ers, such as active Ras, PIP3, Arp2/3, and freshly polymerized actinthat are also concentrated in the inner part of actin waves. On theother hand, markers of the cell back, such as PIP2, myosin II, and cor-texillin, correspond to low values of U and are found outside the wavearea [SDGE + b allows to continuouslyshift from bistable to excitable dynamics, both of which are observedin experiments with D. discoideum cells. In Figure 1D-F, the results of numerical simulations of this model are displayed. Examples forbounded domains that correspond to normal-sized and giant cells areshown, as well as results with periodic boundary conditions that wereused in the subsequent analysis.Model parameters, such as the diffusivities, are typically chosen inan ad hoc fashion to match the speed of intracellular waves with theexperimental observations. The approach introduced in Section 2 nowallows us to estimate diffusivities from data in a more rigorous manner.On the one hand, we may test the validity of our method on in silico data of model simulations, where all parameters are predefined. On theother hand, we can apply our method to experimental data, such as therecordings of cortical actin waves displayed in Figure 1. This will yieldan estimate of the diffusivity of the activator U , as dense areas of fil-amentous actin reflect increased concentrations of activatory signalingcomponents. Note, however, that the estimated value of D U should notbe confused with the molecular diffusivity of a specific signaling mole-cule. It rather reflects an effective value that includes the diffusivitiesof many activatory species of the signaling network and is furthermoreaffected by the specific two-dimensional setting of the model that nei-ther includes the kinetics of membrane attachment/detachment nor thethree-dimensional cytosolic volume.4. Diffusivity Estimation onSimulated and Real Data
In this section we apply the methods from Section 2 to syntheticdata obtained from a numerical simulation and to cell data stemmingfrom experiments as described in Section 3. We follow the formalismfrom Theorem 4 and perform a Fourier decomposition on each dataset. If we set φ k ( x ) = cos(2 πkx ) for k ≤ φ k ( x ) = sin(2 πkx ) for k >
0, then Φ k,l ( x, y ) = φ k ( x/L ) φ l ( y/L ), k, l ∈ Z , form an eigenbasisfor − ∆ on the rectangular domain D = [0 , L ] × [0 , L ]. The corre-sponding eigenvalues are given by λ k,l = 4 π (( k/L ) + ( l/L ) ). Asbefore, we choose an ordering (( k N , l N )) N ∈ N of the eigenvalues (exclud-ing λ , = 0) such that λ N = λ k N ,l N is increasing.In the sequel, we will use different versions of ˆ θ N which correspondto different model assumptions on the reaction term F , concerningboth the effects included in the model and a priori knowledge on theparametrization. While all of these estimators enjoy the same asymp-totic properties as N → ∞ , it is reasonable to expect that they exhibithuge qualitative differences for fixed N ∈ N , depending on how muchknowledge on the generating dynamics is incorporated. In order to de-scribe the model nonlinearities that we presume, we use the notation F , F , F as in (10), (11), (12). As a first simplification, we substitute IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 17 F by (cid:101) F given by (cid:101) F ( U ) = − U ( u − U )(46)in all estimators below. This corresponds to an approximation of thefunction a by an effective average value ¯ a >
0. While this clearly doesnot match full model, we will see that it does not pose a severe restric-tion as a ( U ) tends to stabilize in the simulation. Recall the explicitrepresentation (21) of ˆ θ N . As before, K is the number of nuisance pa-rameters appearing in the nonlinear term F . We construct the followingestimators which capture qualitatively different model assumptions:(1) The linear estimator ˆ θ lin ,N results from presuming K = 0 and F = 0.(2) The polynomial or Schl¨ogl estimator ˆ θ pol ,N , where K = 0 and F ( u ) = k u ( u − u )( u − ¯ au )= − k ¯ au u ( u − u ) + k u ( u − u )= θ (cid:101) F ( u ) + θ F ( u )for known constants k , u , ¯ a > θ = k u ¯ a , θ = k . Thecorresponding SPDE (5) is called stochastic Nagumo equationor stochastic Schl¨ogl equation and arises as the limiting case (cid:15) → full or FitzHugh-Nagumo estimator ˆ θ full ,N , where K = 0and F ( u ) = k u ( u − u )( u − ¯ u ) − k v = θ (cid:101) F ( u ) + θ F ( u ) + θ F ( u )(47) with θ = k u ¯ a , θ = k and θ = k (cid:15)b , where v is given by v t = (cid:82) t e ( t − r )( D V ∆ − (cid:15)I ) u r d r . As before, k , k , u , ¯ a, D v , (cid:15) > θ full ,N in order to estimate different subsets ofmodel parameters at the same time. We use the notation ˆ θ i,N , where i is the number of simultaneously estimated parameters. More precisely,we set ˆ θ ,N = ˆ θ full ,N , and additionally:(1) The estimator ˆ θ ,N results from K = 1 and F θ given by (47)for known θ , θ >
0. This corresponds to an unknown ¯ a .(2) The estimator ˆ θ ,N results from K = 2 and F θ ,θ given by (47).Only θ is known.(3) The estimator ˆ θ ,N results from K = 3 and F θ ,θ ,θ given by(47). All three parameters θ , θ , θ are unknown.In all estimators in this section we set α = 0. This is a reasonablechoice if the driving noise in (7), (8) is close to white noise.It is worth to note that ˆ θ lin ,N is invariant under rescaling the intensityof the data, i.e. substituting X by cX , c >
0. This has the advantage m / s ̂θ lin̂N0 ̂θ pol̂N0 ̂θ full̂N0 m / s ̂θ full̂N0 ̂a=0.1̂θ full̂N0 ̂a=0.2̂θ full̂N0 ̂a=0.30 100 200 300 400 500 600 700 800N−4−3−2−101234 m / s ̂θ ̂θ ̂θ m / s ̂θ ,̂full̂domain̂θ ,̂section̂θ ,̂periodified Figure 2.
Performance of diffusivity estimators on sim-ulated data under different model assumptions in thespatial high-frequency regime. Solid black line is plottedat the true parameter θ = 1 × − , dashed black lineis plotted at zero. In all displays, we restrict to N ≥ Performance on Synthetic Data.
First, we study the perfor-mance of the mentioned estimators on simulations. The numericalscheme is specified in Appendix B. While we have perfect knowledgeon the dynamical system which generates the data in this setting, itis revelatory to compare the different versions of ˆ θ N which correspondto varying levels of model misspecification. The simulation shows that a ( | U t | L ) stabilizes at a value slightly larger than 0 .
1. We demonstratethe effect of qualitatively different model assumptions on our methodin Figure 2 (top left) by comparing the performance of ˆ θ lin ,N , ˆ θ pol ,N and ˆ θ full ,N . The result can be interpreted as follows: As ˆ θ lin ,N doesnot see any information on the wave fronts, the steep gradient at thetransition phase leads to a low diffusivity estimate. On the other hand,ˆ θ pol ,N incorporates knowledge on the wave fronts as they appear in theSchl¨ogl model, but the decay in concentration due to the presence ofthe inhibitor is mistaken as additional diffusion. Finally, ˆ θ full ,N con-tains sufficient information on the dynamics to give a precise estimate. IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 19
In Figure 2 (top right), we show the effect of wrong a priori assump-tions on ¯ a in ˆ θ full ,N . Even for N = 800, the precision of ˆ θ full ,N clearlydepends on the choice of ¯ a . Remember that there is no true ¯ a in theunderlying model, rather, ¯ a serves as an approximation for a ( | U t | L ).Better results can be achieved with ˆ θ ,N , ˆ θ ,N and ˆ θ ,N , see Figure 2(bottom left): ˆ θ ,N has no knowledge on ¯ a and recovers the diffusivityprecisely, and even ˆ θ ,N performs better than the misspecified ˆ θ full ,N from the top right panel of Figure 2.4.2. Discussion of the periodic boundary.
In Figure 2 (bottomright) we sketch how the assumption of periodic boundary conditionsinfluence influences the estimate. While ˆ θ ,N works very well on thefull domain of 200 ×
200 pixels with periodic boundary conditions, itdecays rapidly if we just use a square section of 75 ×
75 pixels. Infact, the boundary conditions are not satisfied on that square section.This leads to the presence of discontinuities at the boundary. Thesediscontinuities, if interpreted as steep gradients, lower the observed dif-fusivity. Hence, a first guess to improve the quality is to mirror thesquare section along each axis and glue the results together. In thismanner we construct a domain with 150 ×
150 pixels, on which ˆ θ ,N performs well. We emphasize that while this periodification proce-dure is a natural approach, its performance will depend on the specificsituation, because the dynamics at the transition edges will still notobey the true underlying dynamics. Furthermore, by modifying thedata set as explained, we change its resolution, and consequently, adifferent amount of spectral information may be included into ˆ θ ,N forinterpretable results.4.3. Performance on Real Data.
A description of the experimen-tal setup can be found in Appendix B. In order to apply the esti-mators from the previous section, further preprocessing is needed, asthe concentration is given in grey values ranging from 0 to 255. Be-fore performing the estimation procedure, we standardize this range to[0 , u ], u = 1, in order to match the stable fixed points of the bistablepolynomial f . Note that this is necessary for all estimators exceptˆ θ lin ,N . In Figure 3 (top left) the behaviour of four of the estimatorson a sample cell is shown. Interestingly, the model-free linear estima-tor ˆ θ lin ,N is close to ˆ θ ,N and ˆ θ ,N , which impose very specific modelassumptions. This pattern can be observed across different cell datasets. In particular, this is notably different from the performance ofthese three estimators on synthetic data. This discrepancy seems toindicate that the lower order reaction terms in the activator-inhibitormodel are not fully consistent with the information contained in theexperimental data. This can have several reasons, for example, it ispossible that a more detailed model reduction of the known signalling m / s lin̂N0 ̂θ ̂θ ̂θ m / s lin̂N0 ̂θ ̂θ ̂θ ,̂T=66.0̂θ ,̂T=132.0̂θ ,̂T=198.0 0 100 200 300 400 500 600 700 800N0.00.20.40.60.81.0 ̂θ ,̂T=258.02̂θ ,̂T=516.04̂θ ,̂T=774.06 Figure 3.
In all displays, we restrict to N ≥
25 inorder to avoid artefacts stemming from low resolution. (top)
Performance of different diffusivity estimators oncell data (top left) and periodified cell data (top right).Dashed black line is plotted at zero. (bottom)
Esti-mation of the reaction parameter θ on simulated data(bottom left) and cell data (bottom right). All times aregiven in seconds.pathway inside the cell is needed. On the contrary, ˆ θ ,N seems to becomparatively rigid in its a priori choices for θ and θ , but eventuallyapproaches the other estimators. Variations in the value of u have animpact on the results for small N but not on the asymptotic behaviour.In Figure 3 (top right), the cell from Figure 3 (top left) is periodifiedbefore evaluating the estimators. As expected from the discussion inSection 4.2, the estimates rise, but the order of magnitude does notchange drastically.4.4. Estimating θ . When solving the linear MLE equations in orderto obtain ˆ θ ,N , we simultaneously get an estimate ˆ θ ,N for θ = k u ¯ a .Note that u = 1 by convention, and θ = k = 1 is treated as knownquantity in this case. Thus, θ is expected to contain information on¯ a . In Figure 3 (bottom) we show the results for ˆ θ ,N for simulateddata (bottom left) and a cell data sample (bottom right). Note that ingeneral, we cannot expect the precision of ˆ θ ,N to increase as N grows,because the reaction term is of order zero. However, it is informativeto consider also the large time regime T → ∞ , i.e. to include moreand more frames to the estimation procedure. In the case of simulateddata, the effective value ¯ a is recovered rather precisely, even for small T ,with increasing precision as T grows. Clearly, this depends heavily onthe model assumptions. In the case of cell data, the results are rather IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 21 m / s lin̂N0 N (X) - inside the cellA N (X) - outside the cell Figure 4.
As before, we restrict to N ≥ (left) Effective diffusivity outside the cell, plot for one dataset. Dashed line is plotted at zero. (right)
Comparisonof the energy inside and outside the cell. Both data setshave the same spatial and temporal extensions.stable. This indicates that it may be reasonable to use the concept of an”effective unstable fixed point ¯ a of the reaction dynamics, conditionedon the model assumptions included in ˆ θ ,N ”, when evaluating cell datastatistically.4.5. The Case of Pure Noise Outside the Cell.
If the data setdoes not contain parts of the cell but rather mere noise, the estima-tion procedure still returns a value. This ”observed diffusivity” (seeFigure 4) originates completely from white measurement noise. Moreprecisely, the appearance and vanishing of singular pixels is interpretedas instantaneous (i.e. within the time between two frames) diffusionto the steady state. Thus, the observed diffusivity in this case can beexpected to be even larger than the diffusivity inside the cell. In thissection, we give a heuristical explanation for the order of magnitude ofthe effective diffusivity outside the cell.We work in dimension d = 2. Assume that a pixel has width x > φ ( y ) with stan-dard deviation σ = x . This way, the inflection points of the (one-dimensional marginal) density match the sharp edges of the pixel. Now,using φ as an initial condition for the heat equation on the wholespace R , the density φ t after time t is also a Gaussian density, withstandard deviation σ t = (cid:112) σ + 2 θt , obtained by convolution with theheat kernel. The maximal value f t max of φ t is attained at y = 0 with f t max = (2 πσ t ) − = (2 π ( σ + 2 θt )) − . Now, if we observe after time t > b >
0, i.e.(48) bf t max ≤ f , this leads to an estimate for the diffusivity of the form(49) θ ≥ ( b − σ t . For example, set t = 0 . s and x = 2 . × − m , as in the data setfrom Figure 4 (left). The intensity decay factor varies between differentpixels in the data set, reasonable values are given for b ≤
30. If b = 30,we get θ ≥ . × − m /s , for b = 20, we get θ ≥ × − m /s ,and for b = 15, we have θ ≥ . × − m /s . The observed diffusivityoutside the cell from Figure 4 is indeed of order 1 × − m s . This givesa heuristical explanation for the larger effective diffusivity outside thecell compared to the estimated values inside the cell.It is important to note that even if the effective diffusivity outsidethe cell is larger, this has almost no effect on the estimation procedureinside the cell. This is because the total energy A N ( X ) , of the noiseoutside the cell is various orders of magnitude smaller than the totalenergy of the signal inside the cell, see Figure 4 (right).5. Discussion and Further Research
In this paper we have extended the mathematical theory of parame-ter estimation of stochastic reaction-diffusion system to the joint esti-mation problem of diffusivity and parametrized reaction terms withinthe variational theory of stochastic partial differential equations. Wehave in particular applied our theory to the estimation of effective dif-fusivity of intracellular actin cytoskeleton dynamics.Traditionally, biochemical signaling pathways were studied in a purelytemporal manner, focusing on the reaction kinetics of the individualcomponents and the sequential order of the pathway, possibly includ-ing feedback loops. Relying on well-established biochemical methods,many of these temporal interaction networks could be characterized.However, with the recent progress in the in vivo expression of fluores-cent probes and the development of advanced live cell imaging tech-niques, the research focus has increasingly shifted to studying the fullspatiotemporal dynamics of signaling processes at the subcellular scale.To complement these experiments with modeling studies, stochasticreaction-diffusion systems are the natural candidate class of modelsthat incorporate the relevant degrees of freedom of intracellular signal-ing processes. Many variants of this reaction-diffusion framework havebeen proposed in an empirical manner to account for the rich plethoraof spatiotemporal signaling patterns that are observed in cells. How-ever, the model parameters in such studies are oftentimes chosen inan ad hoc fashion and tuned based on visual inspection, so that thepatterns produced in model simulations agree with the experimentalobservations. A rigorous framework that allows to estimate the param-eters of stochastic reaction-diffusion systems from experimental datawill provide an indispensable basis to refine existing models, to testhow well they perform, and to eventually establish a new generation of
IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 23 more quantitative mathematical models of intracellular signaling pat-terns.The question of robustness of the parameter estimation problemw.r.t. to specific modelling assumptions of the underlying stochasticevolution equation is an important problem in applications that needsto be further investigated in future research. In particular, this ap-plies to the dependence of diffusivity estimation on the domain and itsboundary. In this work, we based our analysis on a Fourier decompo-sition on a rectangular domain with periodic boundary conditions. Anatural, boundary-free approach is using local estimation techniques asthey have been developed in [AR20] in the context of linear SPDE andgeneralized in [ACP20] to local estimation of the diffusivity in semi-linear stochastic SPDE. An additional approach aiming in the samedirection is the application of a wavelet transform.It is a crucial task to gather further information on the reaction termfrom the data. Principally, this cannot be achieved in a satisfactoryway on a finite time horizon, so the long-time behaviour of maximumlikelihood-based estimators needs to be studied in the context of sto-chastic reaction-diffusion systems. We will address this issue in detailin future work.
Appendix A. Proof of Proposition 1
We prove Proposition 1 by a series of lemmas. First, we show that(7), (8) is well-posed in H = L × L . First note that for any u, v ∈ L ( D ) and x ∈ R : (cid:104) ∂ u f ( x, u ) v, v (cid:105) L (cid:46) | v | L (50)as ∂f ( x, u ) is bounded from above uniformly in x and u . Lemma 9.
There is a unique H -valued solution ( U, V ) to (7) , (8) , and E (cid:20) sup ≤ t ≤ T | ( U t , V t ) | p H (cid:21) < ∞ (51) for any p ≥ . In particular, E (cid:20) sup ≤ t ≤ T | U t | pL (cid:21) < ∞ , E (cid:20) sup ≤ t ≤ T | V t | pL (cid:21) < ∞ (52)for any p ≥ Proof of Lemma 9.
This follows from [LR15, Theorem 5.1.3]. In orderto apply this result, we have to test the conditions ( H , ( H (cid:48) ) , ( H , ( H (cid:48) )(hemicontinuity, local monotonicity, coercivity and growth) therein. Set V = H × H . Define A ( U, V ) = ( A ( U, V ) , A ( U, V )) with A ( U, V ) = D U ∆ U + k f ( | U | L , U ) − k V, (53) A ( U, V ) = D V ∆ V + (cid:15) ( bU − V ) . (54)As B = ( − ∆) − γ is constant and of Hilbert-Schmidt type, we can ne-glect it in the estimates. In order to check the conditions , we haveto look separately at A and A . The statements for A is trivial bylinearity, so we test the parts of the conditions corresponding to A .( H
1) is clear as f ( x, u ) is a polynomial in u and continuous in x . ( H H − (cid:104) A ( U , V ) − A ( U , V ) , U − U (cid:105) H (cid:46) k H − (cid:104) f ( | U | L , U ) − f ( | U | L , U ) , U − U (cid:105) H + k | V − V | L | U − U | L (cid:46) k H − (cid:104) ∂ u f ( | U | L (cid:101) V )( U − U ) , U − U (cid:105) H + k H − (cid:104) ∂ x f (˜ x, U )( | U | L − | U | L ) , U − U (cid:105) H + C (cid:107) ( U , V ) − ( U , V ) (cid:107) H (cid:46) k | ∂ x f (˜ x, U ) | L |(cid:107) U (cid:107) L − (cid:107) U (cid:107) L | | U − U | L + C (cid:107) ( U , V ) − ( U , V ) (cid:107) H (cid:46) (1 + | ∂ x f (˜ x, U ) | L ) (cid:107) ( U , V ) − ( U , V ) (cid:107) H for some ˜ x ∈ R and (cid:101) V : D → R , and | ∂ x f (˜ x, U ) | L = | u a (cid:48) (˜ x ) U ( u − U ) | L (cid:46) | U | L + | U | L . With | U | L = | U | L (cid:46) | U | H , we see thatlocal monotonicity as in ( H
2) is satisfied by taking ρ (( u, v )) = c | u | H ≤ c | ( u, v ) | V . Now, for ( H H − (cid:104) A ( U, V ) , U (cid:105) H ≤ − D U | U | H + k H − (cid:104) f ( | U | L , U ) , U (cid:105) H + k | U | L | V | L (cid:46) − D U | U | H + k (cid:104) ∂ u f ( | U | L , (cid:101) V ) U, U (cid:105) L + C (cid:107) ( U, V ) (cid:107) H (cid:46) − D U | U | H + C | U | L + C (cid:107) ( U, V ) (cid:107) H for some (cid:101) V : D → R , again using (50). Finally, | A ( U ) | H − (cid:46) D U | U | H + k | f ( | U | L , U ) | H − + k | V | L , so it remains to control | f ( | U | L , U ) | H − (cid:46) | f ( | U | L , U ) | L (cid:46) | U | L + | U | L + | U | L , and we have | U | L (cid:46) | U | L as well as | U | L = | U | L and | U | L = | U | L (cid:46) | U | H / (cid:46) | U | H | U | L . Thus ( H
4) is true. Puttingthings together, we get that [LR15] is applicable for sufficiently large p ≥
1, and the claim follows. (cid:3)
In order to improve (52) to the H -norm, it suffices to prove coer-civity, i.e. ( H
3) from [LR15], for the Gelfand triple L ⊂ H ⊂ H . IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 25
Lemma 10.
The solution ( U, V ) to (7) , (8) satisfies ( A ) , i.e. for any p ≥ : E (cid:20) sup ≤ t ≤ T | U t | pH (cid:21) < ∞ , E (cid:20) sup ≤ t ≤ T | V t | pH (cid:21) < ∞ . (55) Proof.
This is done by L (cid:104) A ( U, V ) , U (cid:105) H ≤ − D U | U | H + k (cid:104) ∂ u f ( | U | L , U ) ∇ U, ∇ U (cid:105) L + k | U | H | V | H (cid:46) − D U | U | H + C | U | H + C (cid:107) ( U, V ) (cid:107) H × H , where (50) has been used componentwise. By [LR15, Lemma 5.1.5] weimmediately obtain (55). (cid:3) Remember that by integrating (8), we can write (7) asd U t = ( D U ∆ U t + θ F ( U t ) + θ F ( U t ) + θ F ( U )( t ))d t + B d W t , (56)where we adopted the notation from (10), (11), (12). Lemma 11.
The process ( U, V ) satisfies ( A s ) for some s > .Proof. Let U = U + (cid:101) U be the decomposition of U into its linear andnonlinear part as in (14),(15). We will prove E (cid:20) sup ≤ t ≤ T | (cid:101) U t | pW s,q (cid:21) < ∞ (57)for any p, q ≥ s <
2. As γ > d/ d ≤
2, there is s > U with q = 2, and the claim follows. Let us nowprove (57). First, note that by Lemma 10 and the Sobolev embeddingtheorem in dimension d ≤
2, (57) is true for U t with any p, q ≥ s = 0. For k ∈ N , by | U kt | pL q = | U t | kpL kq , we see that polynomials in U satisfy (57) with s = 0, too. In particular, using that a is bounded, E (cid:20) sup ≤ t ≤ T | F i ( U t ) | pL q (cid:21) < ∞ (58)for i = 1 , p, q ≥
1. A simple calculation shows that the same istrue for F . As in the proof of Proposition 3, we see thatsup ≤ t ≤ T | (cid:101) U t | η,q ≤ | U | η,q + 2 (cid:15) T (cid:15) sup ≤ t ≤ T | θ F ( U ) + θ F ( U ) + θ F ( U ) | L q with η < q ≥ (cid:15) = 2 − η , and the proof is complete. (cid:3) Proof of Proposition 1. In d ≤ H s is an algebra for s >
1, i.e. uv ∈ H s for u, v ∈ H s [AF03]. Together with the assumption that a is bounded, it follows immediately that ( F s,η ) holds for F , F separately(with K = 0) for any s > η <
2. For F , we have for any δ > ≤ t ≤ T | F ( U t ) | s +2 − δ ≤ (cid:90) T | e ( t − r )( D V ∆ − (cid:15)I ) U r | s +2 − δ d r (59) ≤ (cid:90) T ( t − r ) − δ/ | U r | s d r (60) ≤ δ T δ sup ≤ t ≤ T | U t | s , (61)so F satisfies ( F s,η ) even with s ∈ R , η <
4. It is now clear that the( F s,η ) holds for F = θ F + θ F + θ F for any s > η < K = 3). The claim follows from Proposition 3 (i) for U and (ii) for V ,noting that V = (cid:15)bF ( U ). (cid:3) Appendix B. Methods
Numerical simulation.
For the evaluation in Section 4 we simulated(7), (8) on a square with side L = 75, with spatial increment d x =0 . t ∈ [0 , T ], T = 200, with temporal increment d t = 0 . b = 0 . a ( x ) = 0 . − b + 0 . x/ (0 . u L ) − L = 75, with spatial increment d x = 0 .
15, for t ∈ [0 , T ], T = 1000,with temporal increment d t = 0 . =113 µ m and A =2290 µ m . We used b = 0 . a ( x ) = 0 . − b + 0 . x/ (0 . u A ) − D U = 0 . , D V =0 . , k = k = 1 , u = 1 , (cid:15) = 0 . , γ = 0 , σ = 0 .
1. The unit length is1 µ m, the unit time is 1 s . Experimental settings.
For experiments
D. discoideum
AX2 cellsexpressing mRFP-LimE∆ as a marker for F-actin were used. Cell cul-ture, electric pulse induced fusion to create giant cells and live cellimaging were performed as described in [GEW + Acknowledgements
This research has been partially funded by Deutsche Forschungsge-meinschaft (DFG) - SFB1294/1 - 318763901. SA acknowledges fundingfrom MCIU and FEDER – PGC2018-095456-B-I00.
References [ACP20] R. Altmeyer, I. Cialenco, and G Pasemann,
Parameter estimation forsemilinear spdes from local measurements , Preprint: arXiv:2004.14728[math.ST], 2020.
IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 27 [AF03] R. A. Adams and J. J. F. Fournier,
Sobolev Spaces , second ed., Pureand Applied Mathematics, vol. 140, Elsevier/Academic Press, 2003.[AF09] Sarah J. Annesley and Paul R. Fisher,
Dictyostelium discoideumamodel for many reasons , Molecular and Cellular Biochemistry (2009), no. 1-2, 73–91.[AR20] R. Altmeyer and M. Reiß,
Nonparametric Estimation for LinearSPDEs from Local Measurements , Annals of Applied Probability(2020), to appear.[ASB18] Sergio Alonso, Maike Stange, and Carsten Beta,
Modeling randomcrawling, membrane deformation and intracellular polarity of motileamoeboid cells , PLOS ONE (2018), no. 8, e0201977.[ASM +
10] Yoshiyuki Arai, Tatsuo Shibata, Satomi Matsuoka, Masayuki J. Sato,Toshio Yanagida, and Masahiro Ueda,
Self-organization of the phos-phatidylinositol lipids signaling system for random cell migration , Pro-ceedings of the National Academy of Sciences (2010), no. 27,12399–12404.[BAB08] Carsten Beta, Gabriel Amselem, and Eberhard Bodenschatz,
Abistable mechanism for directional sensing , New Journal of Physics (2008), no. 8, 083015.[BBPSP14] Laurent Blanchoin, Rajaa Boujemaa-Paterski, Ccile Sykes, and JuliePlastino, Actin Dynamics, Architecture, and Mechanics in Cell Motil-ity , Physiological Reviews (2014), no. 1, 235–263.[BK17] Carsten Beta and Karsten Kruse, Intracellular Oscillations and Waves ,Annual Review of Condensed Matter Physics (2017), no. 1, 239–264.[BT19] M. Bibinger and M. Trabs, On central limit theorems for power vari-ations of the solution to the stochastic heat equation , Stochastic Mod-els, Statistics and Their Applications (A. Steland, E. Rafaj(cid:32)lowicz, andO. Okhrin, eds.), Springer Proceedings in Mathematics and Statistics,vol. 294, Springer, 2019, pp. 69–84.[BT20] ,
Volatility Estimation for Stochastic PDEs Using High-Frequency Observations , Stochastic Processes and their Applications (2020), no. 5, 3005–3052.[CDVK20] I. Cialenco, F. Delgado-Vences, and H.-J. Kim,
Drift Estimation forDiscretely Sampled SPDEs , Stochastics and Partial Differential Equa-tions: Analysis and Computations (2020), to appear.[CGH11] I. Cialenco and N. Glatt-Holtz,
Parameter estimation for the stochasti-cally perturbed Navier-Stokes equations , Stochastic Process. Appl. (2011), no. 4, 701–724.[CH19] I. Cialenco and Y. Huang,
A Note on Parameter Estimationfor Discretely Sampled SPDEs , Stochastics and Dynamics (2019),doi:10.1142/S0219493720500161, to appear.[Cho19a] C. Chong,
High-Frequency Analysis of Parabolic Stochastic PDEs , An-nals of Statistics (2019), to appear.[Cho19b] ,
High-Frequency Analysis of Parabolic Stochastic PDEs withMultiplicative Noise: Part I , Preprint: arXiv:1908.04145 [math.PR],2019.[Cia18] I. Cialenco,
Statistical inference for SPDEs: an overview , Stat. Infer-ence Stoch. Process. (2018), no. 2, 309–329.[CK20] I. Cialenco and H.-J. Kim, Parameter estimation for discretely sampledstochastic heat equation driven by space-only noise revised , Preprint:arXiv:2003.08920 [math.PR], 2020. [CSS05] John Condeelis, Robert H. Singer, and Jeffrey E. Segall,
THE GREATESCAPE: When Cancer Cells Hijack the Genes for Chemotaxis andMotility , Annual Review of Cell and Developmental Biology (2005),no. 1, 695–718.[DBE +
17] Peter N. Devreotes, Sayak Bhattacharya, Marc Edwards, Pablo A.Iglesias, Thomas Lampert, and Yuchuan Miao,
Excitable Signal Trans-duction Networks in Directed Cell Migration , Annual Review of Celland Developmental Biology (2017), no. 1, 103–125.[DPZ14] G. Da Prato and J. Zabczyk, Stochastic Equations in Infinite Dimen-sions , second ed., Encyclopedia of Mathematics and its Applications,vol. 152, Cambridge University Press, 2014.[FFAB20] Sven Flemming, Francesc Font, Sergio Alonso, and Carsten Beta,
Howcortical waves drive fission of motile cells , Proceedings of the NationalAcademy of Sciences (2020), no. 12, 6330–6338.[FMU19] Seiya Fukushima, Satomi Matsuoka, and Masahiro Ueda,
Excitabledynamics of Ras triggers spontaneous symmetry breaking of PIP3 sig-naling in motile cells , Journal of Cell Science (2019), 12.[GEW +
14] Matthias Gerhardt, Mary Ecke, Michael Walz, Andreas Stengl,Carsten Beta, and Gnther Gerisch,
Actin and PIP3 waves in giantcells reveal the inherent length scale of an excited state , J Cell Sci (2014), no. 20, 4507–4517.[HKR93] M. H¨ubner, R. Khasminskii, and B. L. Rozovskii,
Two Examples ofParameter Estimation for Stochastic Partial Differential Equations ,Stochastic processes (S. Cambanis, J. K. Ghosh, R. L. Karandikar,and P. K. Sen, eds.), Springer, New York, 1993, pp. 149–160.[HLR97] M. Huebner, S. Lototsky, and B. L. Rozovskii,
Asymptotic Properties ofan Approximate Maximum Likelihood Estimator for Stochastic PDEs ,Statistics and Control of Stochastic Processes (Y. M. Kabanov, B. L.Rozovskii, and A. N. Shiryaev, eds.), World Sci. Publ., 1997, pp. 139–155.[HR95] M. Huebner and B. L. Rozovskii,
On asymptotic properties of max-imum likelihood estimators for parabolic stochastic PDE’s , Probab.Theory Related Fields (1995), no. 2, 143–163.[Hue93] M. Huebner,
Parameter Estimation for Stochastic Differential Equa-tions , ProQuest LLC, Ann Arbor, 1993, Thesis (Ph.D.)–University ofSouthern California.[JS03] J. Jacod and A. N. Shiryaev,
Limit Theorems for Stochastic Processes ,second ed., Grundlehren der Mathematischen Wissenschaften, vol. 288,Springer-Verlag, Berlin, Heidelberg, 2003.[KT19] Z. M. Khalil and C. Tudor,
Estimation of the drift parameter forthe fractional stochastic heat equation via power variation , ModernStochastics: Theory and Applications (2019), no. 4, 397–417.[KU19] Y. Kaino and M. Uchida, Parametric estimation for a parabolic lin-ear spde model based on sampled data , Preprint: arXiv:1909.13557[math.ST], 2019.[Lot03] S. Lototsky,
Parameter Estimation for Stochastic Parabolic Equations:Asymptotic Properties of a Two-Dimensional Projection-Based Esti-mator , Stat. Inference Stoch. Process. (2003), no. 1, 65–87.[LR99] S. V. Lototsky and B. L. Rosovskii, Spectral asymptotics of some func-tionals arising in statistical inference for SPDEs , Stochastic Process.Appl. (1999), no. 1, 69–94. IFFUSIVITY ESTIMATION FOR ACTIVATOR-INHIBITOR MODELS 29 [LR00] S. Lototsky and B. L. Rozovskii,
Parameter Estimation for Stochas-tic Evolution Equations with Non-Commuting Operators , Skorokhod’sideas in probability theory (V. Korolyuk, N. Portenko, and H. Syta,eds.), Institute of Mathematics of the National Academy of Science ofUkraine, Kiev, 2000, pp. 271–280.[LR15] W. Liu and M. R¨ockner,
Stochastic Partial Differential Equations: AnIntroduction , Universitext, Springer, 2015.[LS77] R. S. Liptser and A. N. Shiryayev,
Statistics of Random ProcessesI (General Theory) , Applications of Mathematics, vol. 5, Springer-Verlag, New York, 1977.[LS89] R. Sh. Liptser and A. N. Shiryayev,
Theory of Martingales , Mathe-matics and its Applications (Soviet Series), vol. 49, Kluwer AcademicPublishers Group, Dordrecht, 1989.[PS20] G. Pasemann and W. Stannat,
Drift estimation for stochastic reaction-diffusion systems , Electron. J. Statist. (2020), no. 1, 547–579.[PT07] Jan Posp´ıˇsil and Roger Tribe, Parameter Estimates and Exact Vari-ations for Stochastic Heat Equations Driven by Space-Time WhiteNoise , Stoch. Anal. Appl. (2007), no. 3, 593–611.[SDGE +
09] Britta Schroth-Diez, Silke Gerwig, Mary Ecke, Reiner Hegerl, StefanDiez, and Gnther Gerisch,
Propagating waves separate two states ofactin organization in living cells , HFSP Journal (2009), no. 6, 412–427.[Shu01] M. A. Shubin, Pseudodifferential Operators and Spectral Theory , sec-ond ed., Springer-Verlag, Berlin, Heidelberg, 2001.[VWB +
16] Douwe M. Veltman, Thomas D. Williams, Gareth Bloomfield, Bi-Chang Chen, Eric Betzig, Robert H. Insall, and Robert R. Kay,
Aplasma membrane template for macropinocytic cups , Elife (2016),e20085.[Wey11] H. Weyl, ¨Uber die asymptotische Verteilung der Eigenwerte ,Nachrichten von der Gesellschaft der Wissenschaften zu G¨ottingen,Mathematisch-Physikalische Klasse (1911), 110–117. E-mail address ::