# Modeling Diffusion in the Calcium Release Unit

This tutorial we setup and solve for the diffusion of molecules in the electron microscopy derived calcium release unit geometry generated in the [BlendGAMer tutorial](../tutorials/blendertutorial.html).
In order to run this notebook you will need the following packages:
   
- Python environment
- dolfin
- tqdm


```python
from tqdm import tqdm      # Utility to generate progress bars
import dolfin as d         # import dolfin

# Suppress some of the text that FFC dumps to the notebook
import logging
logging.getLogger('FFC').setLevel(logging.WARNING)
```

An easy way to get everything configure a suitable environment is using [conda](https://docs.continuum.io/anaconda/install/).
Using `conda`, the following steps should configure everything you need. 

    conda create -n fenicsproject -c conda-forge fenics
    conda activate fenicsproject
    pip install jupyter tqdm

If you have not completed the [BlendGAMer tutorial](../tutorials/blendertutorial.html), the tetrahedral mesh of the CRU can be [downloaded directly](https://github.com/ctlee/gamer_tutorials/blob/master/data/CRUTetmesh.xml).

## Model of Diffusion With Decay (Strong Formulation)

We will model the dynamics of a molecule with concentration $u$.
The dynamics of $u$ are modeled as

\begin{equation*}
\frac{\partial u}{\partial t} = D\nabla^2 u - \frac{u}{\tau} \enspace \text{in} \enspace \Omega,
\end{equation*}

where $D$ is the diffusion constant, $\tau$ is a decay constant, $\Omega$ is the cytosolic domain, and $t$ is time.

We define boundary conditions

\begin{align}
D(\mathbf{n} \cdot \nabla u) &= J_\text{in} \enspace\text{on}\enspace \partial\Omega_\text{t-tubule}, \\
D(\mathbf{n} \cdot \nabla u) &= 0 \enspace\text{on}\enspace \partial\Omega_\text{other}, \\
\partial\Omega &= \partial\Omega_\text{t-tubule} \cup \partial\Omega_\text{other},
\end{align}

where $J_\text{in}$ is the inward flux on the t-tubule membrane ($\partial \Omega_\text{t-tubule}$).
No flux boundary conditions are applied to all other boundaries ($\partial \Omega_\text{other}$).

Finally we set an initial condition that there is no $u$ in the system at time $t=0$,

\begin{equation*}
u(\mathbf{x},t=0) = 0.
\end{equation*}

We define the values for several physical parameters, initial condition, simulation length, and timestep next.


```python
# Model parameters:
## Physical parameters
D    = 30              # Diffusion coefficient [nm^2/us]
jin  = 0.5             # Inward flux rate [uM*nm/us]
tau  = 500             # Decay timescale [us]

## Initial condition
u0   = d.Constant(0.0) # Initial concentration of u at t=0 [uM]

## Simulation time control
T    = 5000.0          # Total simulation length [us]
dt   =  100.0          # time step size [us]      
t    =    0.0          # Initial time
numsteps = round(T/dt)

print(f"Simulation will run {numsteps} steps")
```

    Simulation will run 50 steps


## Model of Diffusion With Decay (Weak/Variational Formulation)

The finite element method uses a variational formulation of the PDE which allows us to search for solutions from a larger function space. Equations input into FEniCS are done so through their variational formulation; we will now show how to convert our problem into a valid FEniCS input.

First we multiply the governing PDE by a test function, $v$, coming from a function space $V$, integrate over the domain. 

\begin{equation*}
\int_\Omega \frac{\partial u}{\partial t}v \,dx = \int_\Omega D\nabla^2uv \,dx - \int_\Omega \frac{u}{\tau}v \,dx
\end{equation*}

We then apply the divergence theorem,

\begin{equation*}
\int D\nabla^2uv \,dx = -\int_\Omega D\nabla u \cdot \nabla v \,dx + \int_{\partial\Omega} D(\mathbf{n}\cdot\nabla u)v \,ds,
\end{equation*}

to obtain the following variational formulation:

\begin{equation*}
\boxed{\int_\Omega \frac{\partial u}{\partial t}v \,dx = -\int_\Omega D\nabla u \cdot \nabla v \,dx + \int_{\partial\Omega} D(\mathbf{n}\cdot\nabla u)v \,ds -\int_\Omega \frac{u}{\tau}v \,dx}
\end{equation*}

We discretize the time-derivative with a backward Euler scheme due to its simplicity and unconditional stability properties,

\begin{equation*}
\frac{\partial u}{\partial t} \approx \frac{u - u_n}{\Delta t},
\end{equation*}

where $u_n$ represents the (known) solution computed at the previous timestep. 
We further simplify the variational formulation by inserting the Neumann boundary conditions and using the shorthand notation where $\langle \cdot , \cdot \rangle$ represents the inner-product over $\Omega$ and $\langle \cdot , \cdot \rangle_{\partial\Omega_\text{in}}$ represents the inner-product on the boundary $\partial\Omega_\text{in}$. 
Terms are separated onto the left or right hand sides by their dependence on $u$:

\begin{equation*}
\boxed{\langle u,v \rangle + D\Delta t \langle \nabla u,\nabla v \rangle + \frac{\Delta t}{\tau}\langle u,v \rangle =  \Delta t \langle J_\text{in},v \rangle_{\partial\Omega_\text{in}} + \langle u_n,v \rangle}
\end{equation*}

Notice that since $u_n$ is a known value it is not dependent on the unknown, $u$, and therefore is placed on the right hand side. Also, since $D(\mathbf{n}\cdot \nabla u)=0$ on the no-flux boundary those terms drop out of the variational formulation.

In the abstract form this is written as,

\begin{align}
a(u,v) &= \langle u,v \rangle  + D\Delta t \langle \nabla u,\nabla v \rangle + \frac{\Delta t}{\tau} \langle u,v \rangle, \\
L(v) &= \Delta t \langle J_\text{in},v \rangle_{\partial\Omega_\text{in}} + \langle u_n,v \rangle, 
\end{align}

where $a(u,v)$ is a bilinear form and $L(v)$ is a linear functional.

## Solving the Linear Problem Using Dolfin/FEniCS


```python
# Import the mesh and construct a linear Lagrange function space over the mesh
dolfin_mesh = d.Mesh('../data/CRUTetmesh.xml')
V = d.FunctionSpace(dolfin_mesh,'P',1)
```


```python
# Define trial and test functions
u = d.TrialFunction(V)
v = d.TestFunction(V)
# Known solution corresponding to previous timestep. Initialize with initial condition
un = d.interpolate(u0, V)
```


```python
# Recall that in the CRU tutorial we applied the following boundary markings
# Cytosol = 10
# Mitochondria = 20
# T-tubule = 30
# SR = 40

# define integration domains
face_dimension = 2 # we marked faces which have a topological dimension of 2 (triangles)
# Pull out a subdomain of all faces
meshBoundary = d.MeshFunction("size_t", dolfin_mesh, face_dimension, dolfin_mesh.domains())

# Setup "measures" to refer to volume and boundary elements
# Create a measure for the volume
dx = d.Measure('dx', domain = dolfin_mesh)
# Create a measure for boundary 
ds = d.Measure('ds', domain = dolfin_mesh, subdomain_data = meshBoundary)
```

Recalling the abstract form:
\begin{align}
a(u,v) &= \left(1+\frac{\Delta t}{\tau}\right) \langle u,v \rangle  + D\Delta t \langle \nabla u,\nabla v \rangle, \\
L(v) &= \Delta t \langle J_\text{in},v \rangle_{\partial\Omega_\text{t-tubule}} + \langle u_n,v \rangle,
\end{align}

note that the no-flux boundary conditions on $\partial\Omega_\text{other}$ nullify those terms.


```python
# define bilinear and linear forms
a = (1+dt/tau)*u*v*dx + D*dt*d.inner(d.grad(u),d.grad(v))*dx
# ds(30) is the surface of the t-tubule
L = dt*jin*v*ds(30) + un*v*dx
```

### Solve and Save the Simulation


```python
filename = 'solution/dolfinOut.pvd'
vtkfile = d.File(filename)

# store the initial condition
u = d.Function(V)
un.rename("u", "solution")
u.rename("u", "solution")
vtkfile << (un,t)

# Main loop
for idx in tqdm(range(numsteps)):
    # step forward in time
    t += dt 
    # Solve the linear problem
    d.solve(a==L, u)
    # Write the solution for the current timestep to file
    vtkfile << (u,t) 
    # Copy solution to un
    un.assign(u)
```

    100%|██████████| 50/50 [00:51<00:00,  1.03s/it]


Solutions at each timestep are saved in the $\texttt{.vtk}$ format (the collection of snapshots is stored in the $\texttt{dolfinOut.pvd}$ file) and can be visualized using software such as Paraview.

bafd7a801eb5d96da16d1b9e781a4ce4b9fc1f29
ipynb
Jupyter Notebook
notebooks/cru_diffusion.ipynb
RangamaniLabUCSD/gamer_tutorials
810f8f727f7e0a18cfea97340f1c1db428c4b3c5
CC-BY-4.0
notebooks/cru_diffusion.ipynb
RangamaniLabUCSD/gamer_tutorials
810f8f727f7e0a18cfea97340f1c1db428c4b3c5
CC-BY-4.0
notebooks/cru_diffusion.ipynb
RangamaniLabUCSD/gamer_tutorials
810f8f727f7e0a18cfea97340f1c1db428c4b3c5
CC-BY-4.0
2019-12-16T20:41:08.000Z
2019-12-16T20:41:08.000Z
Qwen/Qwen-72B
1. YES
2. YES