Reaction diffusion models

MSE works with models based on parabolic partial differential equation of the following form (1):

(1)\[\begin{split}\begin{equation} \left\{ \begin{array}{lcll} \partial _t U(t,x) &=& \Delta D(x,\Theta)U(t,x)+f(\Theta,U(t,x)) &, \text{ pour $x\in \Omega$, $t\in]0,T]$}\\ U(0,x)&=&U_0(x,\Theta)&, \text{ pour $x\in \Omega$}\\ \vec{\nabla} D(x,\Theta)U(t,x).\vec{n}&=&0&, \text{pour $x\in \delta \Omega$, $t\in]0,T]$}\\ \end{array} \right. \end{equation}\end{split}\]

Where $\Omega$ is the study zone of dimension 1 or 2. The bondary condition on $\delta \Omega$ depends on the study and will be adjustied in case of interaction with the border.

Diffusion term

The diffusion process is represented by the term $\Delta D(x,\Theta)U(t,x)$. In MSE, this part of the model is implemented in the class Diffusion process.

Sources terms

The sources terms are contains in $f(\Theta,U(t,x))$. These terms are implemented in the class Source Term Model, which uses the model’s database. It is also possible to defined new ones using a simple textual description. The relations between compartments are also defined in these terms.

Initial state

The initial state of the system $U_0(x,\Theta)$ is defined in the dynamical system class Dynamical System.

Observation Process

At this step it can be useful to see the MSE simulator has the function (2):

(2)\[\begin{split}\begin{equation} \left\{ \begin{array}{lcl} S&:& \mathbb{R}^p \mapsto ((t,\Omega) \mapsto \mathbb{R}^m)\\ && \Theta \to S(\Theta), \text{ The funtion satisfyng (1) for $\Theta$}. \end{array} \right. \end{equation}\end{split}\]

Let $\mathcal{L}$ in (3) be a function that can represent a likelihood, a log-likelihood, a contrat function. This function can describe an observation process. The implementation is in the class observation Processus.

(3)\[\begin{split}\begin{equation} \left\{ \begin{array}{lcl} \mathcal{L}&:& ((t,\Omega) \mapsto \mathbb{R}^m) \mapsto \mathbb{R}\\ && U \to \mathcal{L}(U), \text{ eg: the likelihood of $U$}. \end{array} \right. \end{equation}\end{split}\]

Optimisation, Maximum likelihood Edstimation

MSE provides algorithms (BFGS, system adjoint, MCMC) to estimate the model parameters. It consists in looking the function (4) as an objetive function. An example is the tutorial optim.

(4)\[\begin{split}\begin{equation} \left\{ \begin{array}{lcl} \mathcal{O}&:& (\mathbb{R}^p) \mapsto \mathbb{R}\\ && \Theta \to \mathcal{O}(\Theta) =\mathcal{L}(S(\Theta)), \text{ eg : the likelihood of $\Theta$}. \end{array} \right. \end{equation}\end{split}\]

Because optimisation requires computing the variation of the likelihood, MSE provides algorithms to perform this task.

Computing variations using adjoint method

In MSE, this algorithm is implemented in the class Adjoint system. It consists in computing $\partial_{\theta_i} \mathcal{O}(\Theta)$.

Assumptions and notations:

assumption: the likelihood in (3) can be written as $\mathcal{L}(u)=\int_Q v(u)(x,t)$ with $Q=\Omega \times [0,T]$,
$H$ is the Hilbert functions $Q\to\mathbb{R}$, with the dot product $<u,v>_H=\int_Q u.v$,
$v(u) \in H$,
$v'(u) \in L(H,H)$ is $\partial _u v(u)$ such that $v(u+\delta u)=v(u)+v'(u)(\delta u) + o(\delta u)$,
assumption: $v'(u)$ can be written as $v'(u)=(\delta u)(x,t)=\tilde v'(u)(x,t) . \delta u(x,t) \in \mathbb{R}$ ,
leading to $\mathcal{L}'(u)(\delta u)=<\tilde v'(u),\delta u>_H$.

Let $S$ play the role of a simulator parameterized by $\Theta$:

\[\begin{split}\begin{array}{llll} S :& \mathbb{R}^p \to H\\ &\Theta \to S(\Theta)=u_\Theta \end{array}\end{split}\]

Linearization of $S$ at $\Theta$:

\[\begin{split}\begin{array}{llll} S'(\Theta) :& \mathbb{R}^p \to H\\ &d \to \lim_{e\to 0}\frac{S(\Theta+ed)-S(\Theta)}{e} \end{array}\end{split}\]

Operator $S'(\Theta)^*$:

\[\begin{split}\begin{array}{llll} S'(\Theta)^* :& H \to \mathbb{R}^p\\ & \text{such that } \forall (x,y)\in H\times\mathbb{R}^p, \langle x,S'(\Theta) (y)\rangle_H=\langle S'(\Theta)^*(x),y\rangle_{\mathbb{R}^p} \end{array}\end{split}\]

computing variation using the implicit function theorem

Link between $\nabla \mathcal{O}(\Theta) \in \mathbb{R}^p$ and $\mathcal{O}'(\Theta) \in \mathbb{R}^{p*}$:

\[\begin{split}\begin{array}{llll} \mathcal{O}'(\Theta):&\mathbb{R}^p \to \mathbb{R}\\ &\epsilon\to \langle \nabla \mathcal{O}(\Theta),\epsilon\rangle_{\mathbb{R}^p}\\ \end{array}\end{split}\]

That is, $\mathcal{O}'(\Theta)^*=\nabla \mathcal{O}(\Theta)$. Taking into account function compositions, with $\delta \Theta \in \mathbb{R}^p$:

(5)\[\begin{split}\begin{array}{lllc} \mathcal{O}'(\Theta)(\delta \Theta) &=\mathcal{L}'(S(\Theta))S'(\Theta)(\delta \Theta)\\ &=\langle \tilde v'(S(\Theta)),S'(\Theta)(\delta \Theta)\rangle_H \end{array}\end{split}\]

We then obtain:

(6)\[\nabla \mathcal{O}(\Theta) = S'(\Theta)^* (\tilde v'(S(\Theta)))\]

Adding the state function

Let $F:\mathbb{R}^p\times H \to H$ such that $S(\Theta)=u_\Theta$ satisfies $F(\Theta,u_\Theta)=0\in H$.

If we apply the implicit function theorem:

\[\partial_u F(\Theta,u_\Theta)\delta u +\partial_\Theta F(\Theta,u_\Theta)\delta \Theta=0\]

and

\[S'(\Theta)=\frac{\delta _u}{\delta_\Theta}=-\partial_uF(\Theta,u_\Theta)^{-1}\partial_\Theta F(\Theta,u_\Theta)\]

By substituting into (5):

\[\begin{split}\begin{array}{lllc} \mathcal{O}'(\Theta)(\delta \Theta)&=\langle \tilde v'(u_\Theta),-\partial_u(F(\Theta,u_\Theta)^{-1}\partial_\Theta F(\Theta,u_\Theta)(\delta \Theta)\rangle_H\\ &=\langle \underbrace{\partial_u(F(\Theta,u_\Theta)^{-1*}\tilde v'(u_\Theta)}_{\large- P},-\partial_\Theta F(\Theta,u_\Theta)(\delta \Theta)\rangle_H \end{array}\end{split}\]

That is, $P+\partial_u(F(\Theta,u_\Theta)^{-1*}\tilde v'(u_\Theta)=0$.

Writing the weak formulation of the adjoint equation, for all $w \in H$:

(7)\[\boxed{\langle \tilde v'(u_\Theta)+\partial _u F(\Theta,u_\Theta)^*p,w \rangle_H=0}\]

$P$ is called the adjoint state.

\[\nabla \mathcal{O}(\Theta)_i=\nabla \mathcal{O}(\Theta).e_i=\langle P,\partial_\Theta F(\Theta,u_\Theta) e_i \rangle_H\]

Details:

$\partial_uF (\Theta, u_\Theta)\in L(H,H)$
$\partial_uF (\Theta, u_\Theta)^*\in L(H,H)$
Such that $\langle \partial_uF (\Theta, u_\Theta)(u),w\rangle_H=\langle u,\partial_uF (\Theta, u_\Theta)^*(w)\rangle_H$

computing variation using the Lagrangian

Let $L$ be:

\[L(\Theta,u,p)=\mathcal{L}(u)+\langle p,F(\Theta,u) \rangle_H\]

If $u=u_\Theta$, then:

\[L(\Theta,u_\Theta,p)=\mathcal{O}(\Theta)\]

Taking the derivative:

\[\mathcal{O}'(\Theta)\delta \Theta=\partial_\Theta L(\Theta,u_\Theta,p)\delta \Theta+\partial_u L(\Theta,u_\Theta,p)\partial_\Theta u_\Theta\delta \Theta\]

We then choose $p$ such that for all $\delta u$:

\[\partial_u L(\Theta,u_\Theta,p)\delta u=0\]

which means:

\[\mathcal{L}'(u_\Theta)\delta u+\langle p,\partial_uF(\Theta,u_\Theta) \delta u \rangle_H=0\]

\[\langle \tilde v '(u_\Theta)+\partial_uF(\Theta,u_\Theta)^*p,\delta u \rangle_H=0\]

which corresponds to the adjoint equation. We then obtain:

\[\mathcal{O}'(\Theta) \delta \Theta=\partial_\Theta L(\Theta,u_\Theta,p)\delta \Theta\]

(8)\[\boxed{\mathcal{O}'(\Theta) \delta \Theta=\langle p,\partial _\Theta F(\Theta,u_\Theta)\delta \Theta \rangle_H}\]

We can then deduce the components of $\nabla \mathcal{O}(\Theta)$ by choosing $\delta \Theta=e_i$.

In the examples, we will explicitly define and solve/simulate system (7) to obtain the vector $\nabla \mathcal{O}(\Theta)$.

Remarks:

The adjoint system (7) does not depend on the model parameterization! This means that it can be computed without considering the parameterization of the model (an important detail for implementation).

Moreover, the left-hand side of the adjoint system, the matrix in equation (7), does not depend on the observation process. This means that the construction of the linear system for simulating the adjoint state depends only on the model.

Equation (8) corresponds to computing an integral over $Q$, regardless of the nature of the observation process, even if observations are discrete in time and/or space.

The computation of $\nabla \mathcal{O}(\Theta)$ reduces to simulating a linear PDE and computing $p$ integrals over $Q$.

Computing variations using the linear tangent equation

comming soon.

References

Interactions between 2 geometrical domains

Exchanges between geometrically different domains, $D1$ and $D2$, can be modeled by adding source terms or fluxes between $D1$ and $D2$. Without interaction, the variational formulation leads to the linear system $M \left( \begin{array}{c} \alpha \\ \beta \end{array} \right)=B$ , $M$ is the matrix:

\[\begin{split}M= \left( \begin{array}{cc} M1 & M1\_2 \\ M2\_1 & M2 \end{array} \right),\end{split}\]

where the first bloc row is about the dynamic on $D1$, the extract diagonal could be zero.

An interaction consists in adding $addM$ to M:

\[\begin{split}addM= \left( \begin{array}{cc} addM1 & addM1\_2 \\ addM2\_1 & addM2 \end{array} \right)\end{split}\]

The following describes various interactions and the corresponding implementation adopted in MSE.

Interaction 2d1d

We first proposed a resolution model with interaction 2d1d, based on the following study: The influence of a line with fast diffusion on Fisher-KPP propagation, https://arxiv.org/abs/1210.3721 .

Let:

$\delta$ a 1d segment,
$\Omega$ a rectangle based on $\delta$,
$V$ is defined in $\Omega$,
$U$ is defined on $\delta$.

The system coupling $\Omega$ and delta is:

(9)\[\begin{split}\begin{equation} \left\{ \begin{array}{lcl} \partial_t U -D\partial_{xx}U &=& \boldsymbol {V}(x,0,t)-U&,\forall t>0,\forall x\in \delta, \\ \partial_t\boldsymbol {V}-d\Delta \boldsymbol {V}&=&\boldsymbol {V}(1-\boldsymbol {V})&,\forall t>0,\forall (x,y)\in \Omega\\ -d\partial_y\boldsymbol {V}(x,0,t)&=&U(x,t)-\boldsymbol {V}(x,0,t) &,\forall t>0,\forall x\in \delta \end{array}\right. \end{equation}\end{split}\]

The interaction consists in adding the 1d sources terme $\boldsymbol {V}(x,0,t)-U(x,t)$ and the 2d flux $U(x,t)-\boldsymbol {V}(x,0,t)$.

Matrices descritpion of the 2d1d interaction

For notation, we assume that 2d domain is $D1$, corresponding to the first block row $(M1, M1\_2)$.Let

\[\phi_i \quad i\in[0,n] \text{ a } \delta \text{ based function} ,\]

\[\psi_i \quad i\in[0,n] \text{ a } \Omega \text{ based function}.\]

The map $I_{2,1}$ makes the link beetween $D2$ and $D1$ dofs : $I_{2,1}(i)$ is the 2d index geometrialy-equals the 1d index i

In the system (9), we write the Green’s formula on adding fluxes in $\Omega$:

\[-\int_\Omega d\Delta V \psi_i = \int_\Omega d\nabla V \nabla \psi_i - \int_\delta d\partial_n V \psi_i\]

Using the 2d1d transfer:

\[\int_\delta d\partial_n V \psi_i += \sum\beta_j\int_\delta\psi_j\psi_i - \sum\alpha_k\int_\delta\phi_k\psi_i\]

This leads to:

\[(addM1)_{i,j} = \int_\delta\psi_j \psi_i \qquad (add1\_2)_{i,k} = -\int_\delta\phi_k \psi_i\]

Note that many $\int_\delta\psi_j \psi_i$ are null. Rewriting it using $mass_{1d}$, we get:

\[\begin{split}\boxed{\begin{equation} \begin{array}{l} (addM1)_{I_{2,1}(i),I_{2,1}(j)} = \int_\delta \psi_{I_{2,1}(j)}\psi_{I_{2,1}(i)} = \int_\delta\phi_j\phi_i = (mass_{1d})_{i,j}\\ (addM1\_2)_{I_{2,1}(i),k} = -\int_\delta\phi_i\phi_k = -(mass_{1d})_{i,k} \end{array} \end{equation}}\end{split}\]

In the system (9), we now Write the variational formulation of the source term in $\delta$:

\[\int_\delta V \phi_i - \int_\delta U \phi_i = \sum \beta_j \int_\delta \psi_j \phi_i - mass_{1d} \alpha\]

If $\psi$ equals $\phi$ on $\delta$, we get:

\[\begin{split}\boxed{\begin{equation} \begin{array}{l} addM2 = mass_{\delta}=mass_{2}\\ (addM2\_1)_{i,I_{2,1}(l)} = -(mass_{\delta})_{i,l}= -(mass_{2})_{i,l} \end{array} \end{equation}}\end{split}\]

Multi-edge connection case

Let $\delta_k$ for $k \in [1,n]$ be edges sharing a vertex at $x_v$. Let $i_k$ be the index in the basis of $\delta_k$ such that $\phi_{i_k}(x_v) = 1$.

The boundary condition at the junction is:

\[-d\partial_x u_k(t,x_v) = u_k(t,x_v) - \sum_{l \ne k} u_l(t,x_v)\]

Matrix additions for multiple edges that must be multiply by $\gamma$:

\[\begin{split} \boxed{\begin{equation} \begin{array}{l} (addMk)_{i_k,i_k} = -(n-1), \forall k \in [1,n] \\ (addMkl)_{i_k,i_l} = 1, \forall (k,l) \in [1,n]^2,\ k \ne l \end{array} \end{equation}}\end{split}\]

This interaction is not implemented because it is equivalent to applying the 2-edge interaction $(n-1)$ times.

Interaction 1d1d

Let $\delta_1$ and $\delta_2$ sharing the same geometry. This interaction consist in allowing the population transfert between them. $U_1$ and $U_2$ are respecrtively defined on $\delta_1$ and $\delta_2$. The intercation in adding the terms:

(10)\[\begin{split}\begin{equation} \left\{ \begin{array}{lcll} \partial _t U_1(t,x) &+=& \lambda(U_2(t,x)-U_1(t,x)) &, \text{ pour $x\in \delta_1 $, $t\in]0,T]$}\\ \partial _t U_2(t,x) &+=& \lambda(U_1(t,x)-U_2(t,x)) &, \text{ pour $x\in \delta_2 $, $t\in]0,T]$}\\ \end{array} \right. \end{equation}\end{split}\]

Matrices description of the 1d1d interaction

Let

\[\phi_i \quad i\in[0,n] \text{ a } \delta_1 \text{ based function} ,\]

\[\psi_i \quad i\in[0,n] \text{ a } \delta_2 \text{ based function}.\]

Let $I_{\delta1\delta2}(i)$ be the index in $\delta_2$ geometrically equal to the $\delta_1$ index $i$, $\phi_i = \psi_{I_{\delta1\delta2}(i)}$ , $I_{\delta2\delta1}(I_{\delta1\delta2}(i)) = i$.

The weak formulation in $\delta_1$ of the source terms described in (10), leads to :

\[\int_{\delta_1} U_2 \cdot \phi_i - \int_{\delta_1} U_1 \cdot \phi_i = \sum \beta_j \int_{\delta_1} \psi_j \cdot \phi_i - \sum \alpha_j \int_{\delta_1} \phi_j \cdot \phi_i\]

$addM1_{i,j} = \int_{\delta_1} \phi_j \cdot \phi_i = (mass_1)_{i,j}$
$addM1M2_{i,j} = \int_{\delta_1} \psi_j \cdot \phi_i = \int_{\delta_1} \phi_{I_{\delta2\delta1}(j)} \cdot \phi_i = (mass_1)_{i, I_{\delta2\delta1(j)}}$

Thus:

\[\begin{split}\boxed{\begin{equation} \begin{array}{l} addM1_{i,j} = (mass_1)_{i,j}\\ (addM12)_{i,I_{\delta1\delta2(j)}} = -(mass_1)_{i,j} \end{array} \end{equation}}\end{split}\]

In $\delta_2$, we get:

\[\begin{split}\boxed{\begin{equation} \begin{array}{l} addM2_{i,j} = (mass_2)_{i,j}\\ (addM21)_{i,I_{\delta2\delta1(j)}} = -(mass_2)_{i,j} \end{array} \end{equation}}\end{split}\]

Interaction 2d2d

The connexion between two field $(\Omega_1,\Omega_2)$ sharing an border $\delta$ is defined as follow:

\[\begin{split}\begin{equation} \left\{ \begin{array}{lcll} \vec{\nabla} D_1(x,\Theta)U_1(t,x).\vec{n}&+=&\gamma(U_2(x,t)-U_1(x,t))&, \text{pour $x\in \delta $, $t\in]0,T]$}\\ \vec{\nabla} D_2(x,\Theta)U_2(t,x).\vec{n}&+=&\gamma(U_1(x,t)-U_2(x,t))&, \text{pour $x\in \delta $, $t\in]0,T]$}\\ \end{array} \right. \end{equation}\end{split}\]

The tutorial of interaction 2d2d model is avaible here: Tutorial: simulation with 2D2D interaction.

Matrices description of the 2d2d interaction

hese interations are implemented in the class Interactions.