\documentclass[10pt]{article}
\usepackage{a4wide}
\usepackage{listings}
\usepackage{graphicx}
\usepackage{times}
\usepackage{amsfonts}
\usepackage{amssymb}
\newcommand{\iint}[1]{\int\!\!\!\!\int\limits_{\!\!\!#1}}
\newcommand{\iiint}[1]{\int\!\!\!\!\int\!\!\!\!\int\limits_{\!\!\!\!\!\!\!\!\!#1}}
\newcommand{\text}[1]{\mathrm{#1}}

\begin{document}

\section*{Format for delivery of report and programs}
%
The format of the project is that of a printed file or hand-written report.
The programs should also be included with the report. Write only your candidate number on the first
page of the report and state clearly that this is your report for project 5 of FYS3150, fall 2011. 
There will be a box marked 'FYS3150' at the reception of the Department of Physics (room FV128). 

\section*{Project 5, Diffusion of neurotransmitters in the synaptic cleft, deadline  December 12, 3pm}
The dominant way of transporting signals between neurons (nerve cells)
in the brain is by means of diffusion of particular signal molecules
called \emph{neurotransmitters} across the synaptic cleft separating
the cell membranes of the two cells. A drawing of a synapse is
given in Fig.~\ref{fig:figure1}.
\begin{figure}[thb]
\centerline{\includegraphics[width=9cm]{thompsonB2000-p38.eps}}
\caption{\small Drawing of a synapse. The axon terminal is the knoblike
structure and the spine of the receiving neuron is the bottom one. The
synaptic cleft is the small space between the presynaptic (axon)
and postsynaptic (dendritic spine) membrane.
(From Thompson: ``The Brain'', Worth Publ., 2000)}
\label{fig:figure1}
\end{figure}

Following the arrival of an action potential in the axon terminal a
process is initiated in which (i) vesicles inside the axon terminal
(filled with neurotransmitter molecules) merge with the presynaptic
(axon) membrane and (ii) release neurotransmitters into the synaptic
cleft. These neurotransmitters diffuse across the synaptic cleft to
receptors on the postsynaptic side which ``receives'' the signal.
A schematic illustration of this process is shown in
Fig.~\ref{fig:figure2}(left).
\begin{figure}[t]
\centerline{\includegraphics[width=10cm]{thompsonB2000-p39.eps}
\includegraphics[width=3cm]{kandel-B1991-p217-glutamate.eps}
\includegraphics[width=3cm]{kandel-B1991-p217-GABA.eps}}
\caption{\small Left: Schematic drawing
of the process of vesicle release from the axon terminal and release of
transmitter molecules into the synaptic cleft. (From Thompson: ``The
Brain'', Worth Publ., 2000). Right: Molecular structure of the
two important neurotransmitters \emph{glutamate} and \emph{GABA}.}
\label{fig:figure2}
\end{figure}
Since the transport process in the synaptic cleft is governed by
diffusion, we can describe it mathematically by
%
\begin{equation}
\frac{\partial u}{\partial t} = D \nabla^2 u,
\label{eq:diffusion_eq_3D}
\end{equation}
%
where $u\,$ is the concentration of the particular neurotransmitter, and
$D$ is the diffusion coefficient of the neurotransmitter in this
particular environment (solvent in synaptic cleft).

If we assume (i) that the neurotransmitter is released
roughly equally on the ``presynaptic'' side of the synaptic cleft, and
(ii) that the synaptic cleft is roughly equally wide across the whole
synaptic terminal, we can, given the large area of the synaptic cleft
compared to its width, assume that the neurotransmitter concentration
only varies in the direction across the synaptic cleft (from
presynaptic to postsynaptic side). We choose this direction to be the
$x$-direction (see Fig.~\ref{fig:figure3}).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FIGURE : synaptic cleft
\begin{figure}[b]
\centerline{\includegraphics[width=10cm]{synaptic_cleft.eps}}
\caption{\small Schematic drawing of the synaptic cleft in our model. The
black dots represent neurotransmitter molecules, and the situation
shown corresponds to the situation immediately after neurotransmitter
release into the synaptic cleft.}
\label{fig:figure3}
\end{figure}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
In this case $u({\bf r})=u(x)$, the diffusion equation
reduces to
%
\begin{equation}
\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}.
\label{eq:diffusion_eq_1D}
\end{equation}\newline
Immediately after the release of a neurotransmitter into the
synaptic cleft ($t=0$) the concentration profile in the $x$-direction
is given by
%
\begin{equation}
u(x,t=0) = N \, \delta(x),
\label{eq:initial_condition}
\end{equation}
%
where $N$ is the number of particle released into the synaptic cleft
per area of membrane.

To get an idea over the time-dependence of the neurotransmitter
concentration at the postsynaptic side ($x=d$), we can look at the
solution of a ``free'' random walk (i.e., no obstacles or particle
absorbers in either direction). The
solution of Eq.~(\ref{eq:diffusion_eq_1D}) with the initial condition
in Eq.~(\ref{eq:initial_condition}) is given by 
(see Nelson: \emph{Biological Physics}, p. 143 or Lectures notes chapter 12.3)
%
\begin{equation}
u(x,t) = \frac{N}{\sqrt{4 \pi D t}} e^{-x^2/4Dt}\;\;.
\label{eq:solution_delta_1D}
\end{equation}\newline
The concentration at the postsynaptic side $u(d,t)$
approaches 0 in the limit $t \rightarrow 0\;$ and
$t \rightarrow \infty$. 

The above assumption regarding the
neurotransmitter molecules undergoing a ``free'' random walk, is
obviously a simplification. In the true diffusion process in the
synaptic cleft the neurotransmitter molecules will, for example,
occasionally bump into the presynaptic membrane they came from. Also
at the postsynaptic side the neurotransmitters are absorbed by
receptors located on the postsynaptic cell membrane and are thus
(temporally) removed from the solution.

To approach this situation in our mathematical model we can impose the
following boundary and initial conditions with $x\in[0,d]$
%
\begin{equation}
u(x=0,t>0) = u_0, \;\;u(x=d,\mbox{all $t$})=0,
\;\;u(0 < x < d,t < 0) = 0 \;\;.
\label{eq:initial_conditions_2}
\end{equation}
Hereafter we set $d=1$.
This corresponds to that (i) for $t<0$ there are no neurotransmitters
in the synaptic cleft, (ii) for $t>0$ the concentration of
neurotransmitters at the presynaptic boundary of the synaptic
cleft ($x=0$) is kept
\emph{fixed}  at $u=u_0=1$ in our case, and (iii) that the postsynaptic receptors
immediately absorb nearby neurotransmitters so that $u=0$ on the
postsynaptic side of the cleft ($x=d=1$).


The full solution of the diffusion equation with
boundary/initial conditions in Eq.~(\ref{eq:initial_conditions_2})
can be found in a closed form. We will use this solution to test our numerical calculations.


We are thus looking at a one-dimensional
problem 
\[
 \frac{\partial^2 u(x,t)}{\partial x^2} =\frac{\partial u(x,t)}{\partial t}, t> 0, x\in [0,d]
\]
or 
\[
u_{xx} = u_t,
\]
with initial conditions, i.e., the conditions at $t=0$, 
\[
u(x,0)= 0 \hspace{0.5cm} 0 < x < d
\]
with $d=1$ the length of the $x$-region of interest. The 
boundary conditions are 
\[
u(0,t)= 1 \hspace{0.5cm} t > 0,
\]
and 
\[
u(d,t)= 0 \hspace{0.5cm} t > 0.
\]

In this project we want to study the numerical stability of three methods for partial differential equations
(PDEs). 
These methods are 
\begin{enumerate}
\item The explicit forward Euler algorithm with discretized versions of time given by a forward formula and
a centered difference in space resulting in
 \[
u_t\approx \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}=\frac{u(x_i,t_j+\Delta t)-u(x_i,t_j)}{\Delta t}
\]
and
\[
u_{xx}\approx \frac{u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t)}{\Delta x^2},
\]
or
\[
u_{xx}\approx \frac{u(x_i+\Delta x,t_j)-2u(x_i,t_j)+u(x_i-\Delta x,t_j)}{\Delta x^2}.
\]
\item The implicit Backward Euler with
 \[
u_t\approx \frac{u(x,t)-u(x,t-\Delta t)}{\Delta t}=\frac{u(x_i,t_j)-u(x_i,t_j-\Delta t)}{\Delta t}
\]
and
\[
u_{xx}\approx \frac{u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t)}{\Delta x^2},
\]
or
\[
u_{xx}\approx \frac{u(x_i+\Delta x,t_j)-2u(x_i,t_j)+u(x_i-\Delta x,t_j)}{\Delta x^2},
\]
\item Finally we use the implicit Crank-Nicolson scheme with 
a time-centered scheme at $(x,t+\Delta t/2)$
 \[
u_t\approx \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}=\frac{u(x_i,t_j+\Delta t)-u(x_i,t_j)}{\Delta t}.
\]
The corresponding spatial second-order derivative reads
\[
u_{xx}\approx \frac{1}{2}\left(\frac{u(x_i+\Delta x,t_j)-2u(x_i,t_j)+u(x_i-\Delta x,t_j)}{\Delta x^2}+\right.
\]
\[
\left. \frac{u(x_i+\Delta x,t_j+\Delta t)-2u(x_i,t_j+\Delta t)+u(x_i-\Delta x,t_j+\Delta t)}{\Delta x^2}
\right).
\] 
Note well that we are using a time-centered scheme wih $t+\Delta t/2$ as center.
\end{enumerate}
 \begin{enumerate}
\item[a)] Find the closed form solution to this problem. You will need this in order to study the numerical accuracy of your results.  To find the closed-form solution, we will need the  
stationary solution (steady-state solution). The solution to the steady-state problem is on the form $u(x)=Ax+b$. The solution for the steady-state case $u_s$ that obeys the above boundary conditions is 
\[
u_s(x) = 1-x. 
\]
You can use this solution to  define a new function $v(x)=u(x)-u_s(x)$ with boundary conditions
$v(0)=v(d)=0$. The latter is easier to solve both numerically and on a closed form. 
\item[b)] Write down the algorithms for these three methods and the equations you need to implement.
For the implicit schemes show that the equations lead to a tridiagonal matrix system for the new values.
\item [c)] Find the truncation errors of these three schemes and investigate their stability properties.
\item [d)] Implement the three algorithms in the same code and perform tests of the solution 
for these three approaches
for $\Delta x=1/10$, $h=1/100$ using  $\Delta t$ as dictated by the stability limit of the explicit scheme.
Study the solutions at two time points $t_1$ and $t_2$ where $u(x,t_1)$ is smooth but still significantly curved
and $u(x,t_2)$ is almost linear, close to the stationary state.
Remember that for solving the tridiagonal equations you can use your code from project 1.  
\item[e)] Compare the solutions at $t_1$ and $t_2$ with the analytic result for the continuous problem.
Which of the schemes would you classify as the best?
\item[f)]  This part is optional but gives you an additional 20\% on the final score!  
The above problem can be solved using Monte Carlo methods and random walks. We follow here
Farnell and Gibson in Journal of Computational Physics, volume {\bf 208}, pages 253-265 (2005).
Choose a constant step length $l_0=\sqrt{2D\Delta t}$ and equal probability for jumping left and right.
Set up the algorithm for solving the above diffusion problem and write a code to do it.
Compare your results with those from the partial differential equation solution and comment the results.
\item[g)] This part is also optional, and together with f) gives an additional 10\% score. 
Change the above stepsize by using a Gaussian distribution with mean value $1$ and standard deviation $0$. The step length of the random walker is now $l_0=\sqrt{2D\Delta t}\xi$, where $\xi$ is random number chosen from the above Gaussian distribution.  Implement this stepsize to the program from 
f) and compare the results and comment.  
\end{enumerate}



\end{document}








