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. In order to run this notebook you will need the following packages:

Python environment
dolfin
tqdm

[1]:

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. 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, the tetrahedral mesh of the CRU can be downloaded directly.

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.

[2]:

# 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¶

[3]:

# 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)

[4]:

# 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)

[5]:

# 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.

[6]:

# 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¶

[7]:

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.