\documentclass{kluwer}


\usepackage{graphicx}

%\usepackage{pdfsync}

\newcommand{\ee}{\end{equation}}
\newcommand{\etal}{{\it et al.}\ }
\newcommand{\be}{\begin{equation}}
\newcommand{\ea}{\end{eqnarray}}
\newcommand{\ba}{\begin{eqnarray}}
\newcommand{\R}[1]{$(\ref{#1})$}
\newcommand{\rmd}{{\rm d}}
\newcommand{\rme}{{\rm e}}
\newcommand{\rmi}{{\rm i}}
\newcommand{\e}{{\rm e}}
\newcommand{\QQ}{{\cal Q}}

\def\Zset{\Bbb Z}
\def\Rset{\Bbb R}
\def\Tset{\Bbb T}
\def\Cset{\Bbb C}

\newcommand{\PD}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\FD}[2]{\frac{\rmd #1}{\rmd #2}}


\begin{document}
\begin{article}



\begin{opening}
\title{Spatio-temporal filtering properties of a dendritic cable with
  active spines: a modeling study in the spike-diffuse-spike
  framework} 


\author{Y Timofeeva \email{yulia@ma.hw.ac.uk} and G J Lord \email{gabriel@ma.hw.ac.uk}}
\institute{Department of Mathematics, Heriot-Watt University, Edinburgh, EH14 4AS, UK.}

\author{S Coombes \email{stephen.coombes@nottingham.ac.uk}}
\institute{School of Mathematical Sciences, University of Nottingham,
  Nottingham, NG7 2RD, UK.} 



\begin{abstract}

The spike-diffuse-spike (SDS) model describes a passive dendritic tree
with active dendritic spines. Spine-head dynamics is modeled with a
simple integrate-and-fire process, whilst communication between spines
is mediated by the cable equation. In this paper we develop a
computational framework that allows the study of multiple spiking
events in a network of such spines embedded on a simple
one-dimensional cable. In the first instance this system is shown to
support saltatory waves with the same qualitative features as those
observed in a model with Hodgkin-Huxley kinetics in the spine-head. Moreover, there is excellent agreement with the analytically
calculated speed for a solitary saltatory pulse. 
Upon driving the system with time-varying external input we find
that the distribution of spines can play a crucial role in
determining spatio-temporal filtering properties. In particular, the
SDS model in response to periodic pulse train shows a positive
correlation between spine density and low-pass temporal filtering
that is consistent with the experimental results of Rose and Fortune [1999, `Mechanisms for generating temporal filters in the electrosensory system'. {\em The Journal of Experimental Biology} {\bf 202}, 1281--1289].
Further, we demonstrate the robustness of observed wave properties to
natural sources of noise that arise both in the cable and the
spine-head, and highlight the possibility of purely noise induced
waves and coherent oscillations.


\end{abstract}

\keywords{spike-diffuse-spike, dendritic spines, filtering, noise}

\end{opening}



\section{Introduction}

In 1891 Ram\'on y Cajal showed that dendritic spines are present in the dendrites of many neurons of the cerebral cortex of mammals \cite{Cajal1891}.  Dendritic spines are small mushroom like appendages  with a bulbous head and a tenuous stem (of length around $1 \mu$m) and may be found in their hundreds of thousands on the dendritic tree of a single cortical
pyramidal cell.  These extensions of the dendritic tree provide junction points for the axons of other neurons (i.e., provide surface area for synapses), and thus serve as loci for receiving inputs. In the cerebral cortex approximately 80\% 
of all excitatory synapses are made onto dendritic spines. Since the biophysical properties of spines can be modified by experience in response to patterns of chemical and electrical activity, morphological and electro-chemical changes in populations of dendritic spines are thought to provide a basic mechanism for Hebbian learning in the nervous system. 
In fact as far back as 1899 Cajal (\citeyear{Cajal1899}) was arguing
that spines could be involved in learning and that physical changes in
spines were associated with neuronal function, suggesting that they
might grow with activity and retract during inactivity or sleep. 
These notions continue to be influential in the study of the function
of dendritic spines even today \cite{Yuste01,Yuste01a}. 

In recent years the properties of spines have also been linked with
the implementation of logical computations \cite{Shepherd87}
coincidence detection \cite{Larkum99} orientation tuning in complex
cells of visual cortex \cite{Mel98} and the amplification of distal
synaptic inputs \cite{Miller85}.   
At the organismal level there is now evidence to suggest that the
density of dendritic spines may reflect overall mental agility
\cite{Zito02}.  Conversely, many neurological diseases resulting in
mental retardation have been associated with spine loss or spine
morphology changes (such as Fragile
X-syndrome\footnote{http://www.fraxa.org/}).  However, the focus of
this paper will be on the implication of excitable channels in the spine-head membrane for single neuron dynamics. 
The benefits of excitable membrane for the amplification of excitatory
synaptic inputs was first discussed by Jack et. al. (\citeyear{Jack75}).
If dendritic spines possess excitable membrane, the spread of current from one spine along the dendrites may bring adjacent spines to threshold for impulse generation, resulting in a saltatory propagating wave in the distal dendritic branches \cite{Shepherd85}.  However, it is only relatively recently that confocal and two-photon microscopy observations  have confirmed the generation of action potentials in the dendrites (see \cite{Segev98} for a perspective).  


The first step towards the development of a spiny dendritic tissue
model that might be used to explore these issues can be attributed to
Baer and Rinzel (\citeyear{Baer91}) who considered a passive uniform unbranched
dendritic tree coupled to a population of excitable dendritic spines.
In this \textit{continuum} model the active spine-head dynamics is
modeled with Hodgkin-Huxley (HH) kinetics whilst the (distal) dendritic
tissue is modeled with the cable equation.  The spine-head is coupled
to the cable via a spine-stem resistance  that delivers a current
proportional to the number of spines at the contact point.  There is
no direct coupling between neighboring spines; voltage spread by
diffusion along
the cable is the only way for spines to interact.  Although the
numerical studies of Baer and Rinzel (\citeyear{Baer91}) show traveling wave solutions, the
underlying continuous nature of the model precludes the possibility
that these waves are truly saltatory.  The saltatory nature of a
propagating wave in a spiny neuron may be directly attributed to the
fact that active spines are physically separated. Although we can
numerically simulate  the nonlinear and nonuniform properties of
biologically realistic dendritic trees with discrete and clustered
distributions of spines, based around natural extensions of the
Baer-Rinzel (BR) model, there is a lack of analytical tools for
dealing with such systems. 
However, recent work by Coombes, Bressloff and Lord
\cite{Coombes00,Coombes01,Lord02,Coombes03} has shown that the active membrane
dynamics of spines can be treated using an analytically tractable
integrate-and-fire (IF) process.  The resulting model has been termed
the Spike-Diffuse-Spike (SDS) model since spine-head dynamics is an
all-or-nothing action potential response, whilst the dendritic cable is modeled as 
a passive structure.  Not only can saltatory wave propagation be 
naturally analyzed in the SDS framework, the model is computationally inexpensive
and ideally suited for the study of neural response to complicated
spatio-temporal patterns of synaptic input that typically occur in
cortical neurons.  It is precisely these points we choose to highlight in this paper.

In section \ref{SDS} we review the SDS framework and extend previous work by showing how one may express solutions
in terms of a Dyson-like series expansion.  Moreover, we discuss the
numerical implementation of analytical solutions to the SDS model that
depend upon numerically determined firing events.  For comparison with
more conventional (numerical PDE) approaches we also develop an
implementation of the SDS model within the {\sc neuron} simulation
environment \cite{Hines}. As a first illustration of the usefulness of the SDS
model of spiny dendritic tissue we present a case study of saltatory
propagating waves in a model with regularly spaced spines.  Also in
section \ref{saltatory} we treat wave propagation (and its failure)
for more irregular distributions of spines.  After establishing the
ability of the SDS model to accurately describe the sorts of saltatory
waves one finds in the more biophysically detailed BR model,
we turn next, in section \ref{filtering}, to the issue of the active
dendritic tree as a nonlinear filter.  In particular we use the SDS
model to address the observation of Rose and Fortune (\citeyear{Rose99})
that there is a positive correlation between spine density and
low-pass temporal filtering.  In section \ref{noise} we test the
robustness of observed wave and filtering properties to
natural sources of noise that arise both in the cable and the 
spine-head, and highlight the possibility of purely noise induced
waves and coherent oscillations. Finally, in section
\ref{discussion}, we look forward and describe some of the open
problems in single neuron function that may be readily addressed with a
further analysis of SDS dynamics. 


\section{The spike-diffuse-spike model}
\label{SDS}



We consider a uniform passive dendritic cable with a given
distribution of spines along its length. 
A schematic diagram of the SDS model is shown in Figure \ref{schema}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[th]
\begin{center}
\includegraphics[width=3.5in]{SDSdiagram.eps}
\caption{A schematic representation of the SDS model showing a passive cable with electrically connected active spines.}
\label{schema}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
The dynamics of membrane voltage in the cable $V=V(x,t)$ is described by the equation
\be
\pi a C_m \PD{V}{t}=\frac{\pi a^2}{4 R_a}\PD{^2 V}{x^2}-\frac{\pi a}{R_m}V+\rho(x)I_{\rm{sp}}.
\label{cable}
\ee
Here $a$ is a diameter of the cable (measured in $\mu$m), $R_a$ is the specific cytoplasmic
resistivity (in $\Omega \cdot$cm), $C_m$  and $R_m$  are respectively 
the specific membrane capacity (in $\mu$F/cm$^2$)
and a
resistance across a unit area of passive membrane (in $\Omega \cdot$cm$^2$).
Spines are connected to the cable at the discrete
points $x_n$ with the distribution function 
$\rho(x)=\sum_{n\in  \Gamma}\delta(x-x_n)$, where $\Gamma$ is a
discrete set that indexes the spines. 
Each spine generates a sequence of action potentials in its spine-head given by the function $\widehat V(x,t)$.
As a result the spine that just ``fired'' passes the spine current $I_{\rm{sp}}=(\widehat{V}-V)/r$ into the cable.  
The spine stem resistance of an individual spine is given by $r$.
Denoting
$T_n^m,~m\in \Zset$ as the time of the $m$th firing event of
the $n$th spine the function $\widehat{V}(x_n,t)$ is given by
$\widehat{V}(x_n,t)=\sum_m \eta(t-T_n^m)$. Here $\eta(t)$ specifies the universal shape of an action potential.
The generator of action potentials in the $n$th spine-head
evolves according to\ba
\widehat{C} \PD{U_n}{t}=
-\frac{U_n}{\widehat{r}}+\frac{V_n-U_n}{r}-\underbrace{\widehat{C}h\sum_m\delta(t-T_n^m)}_{\rm {reset}},
\label{sp}
\ea
where $V_n=V(x_n,t)$. Here the parameters $\widehat{C}$ and
$\widehat{r}$ describe the electrical properties of the spine-head
membrane, namely its capacitance and resistance respectively.
The spine's firing times $T_n^m$ are defined 
in terms of the IF process according to
\begin{equation}
T_n^m={\rm inf} \{t~|~ U_n(t)\geq h,~ t> T_n^{m-1}+\tau_R\} .
\end{equation}
Hence,  a spine fires whenever $U_n$, driven by current from the shaft,
crosses some threshold potential $h$.
Just after the firing event the variable $U_n$ resets to zero,
modeled by the last term in equation (\ref{sp}). Multiple spiking
events from an individual spine are controlled by a refractory time-scale $\tau_R$ during which the spine can not fire.

It is convenient to rewrite  (\ref{cable}) in the standard form
\be
\PD{V}{t}=D \PD{^2 V}{x^2}-\frac{V}{\tau}+D r_a\rho(x)\frac{\widehat{V}-V}{r},
\label{cable1}
\ee
where $D=\lambda^2/\tau$ denotes the diffusion coefficient of the cable, $\tau=C_m R_m$ is the membrane time constant, $\lambda=\sqrt{aR_m/4 R_a}$ is the electronic space constant and $r_a=4R_a/\pi a^2$ denotes the intracellular resistance per unit length of cable.
The integration of this equation
gives us an implicit expression for the membrane potential in the form 
\be
V(x,t)=\frac{D r_a}{r}\sum_{k \in \Gamma}\int_0^t {\rm d} s~ G(x-x_k,t-s)[\widehat{V}(x_k,s)-V(x_k,s)].
\label{fullsol}
\ee
Here, $G(x,t)$ is the Green's function of the infinite uniform passive cable 
\be
G(x,t)=\frac{1}{\sqrt{4\pi D t}}~{\rm e}^{-\varepsilon t}~{\rm e}^{-x^2/(4Dt)}\Theta(t) ,
\ee
where
$\varepsilon=1/\tau$ and $\Theta(t)$ is the Heaviside step function.
The expression for the cable voltage (\ref{fullsol}) has a Dyson-like form \cite{Bressloff97} suggesting a Neumann series  solution, obtained by repeated substitution of  (\ref{fullsol}) into itself.  Introducing the parameter $\Lambda=Dr_a/r$, and generating just the first two terms in this expansion gives
\ba
\label{firstorder}
&&V(x,t)=
 \Lambda \sum_{k}\int_0^t {\rm d} s
  G(x-x_k,t-s)\widehat V(x_k,s) \\
&&  - \Lambda^{2} \sum_{k,p}\int_0^t {\rm d}s G(x-x_k,t-s)
 \int_0^s {\rm d} s'
G(x_k-x_p,s-s')\widehat{V}(x_p,s') ,\nonumber
\ea
where $\{p,k\} \in \Gamma$.
Under the approximation that $\Lambda \ll 1$, only the first term
significantly contributes to the full solution.  This first term is an
exact solution of the SDS model in the presence of a partial current
flux $I_{\rm sp}=\widehat V/r$ from the spine-head into the cable
instead of the original full flux. This gives rise to the so-called {\em partial SDS model}.
In this case the solution takes the explicit form 
\be
V(x,t)=\Lambda \sum_{k \in \Gamma,m}H(x-x_k,t-T_k^m),\qquad \max_{k,m}\{T_k^m\}\leq t< T^\ell_j,
\label{sol}
\ee 
where $m=m(k)$ counts firing events at each spine and $H(x,t)=\int_0^t G(x,t-s)\eta(s){\rm d}s$.  
For a simple action potential shape given  
by a rectangular pulse $\eta(t)=\eta_0\Theta(t)\Theta(\tau_S-t)$ (with strength $\eta_0$ and duration $\tau_S$), the function $H(x,t)$ can be found in closed form \cite{Coombes03} 
as $H(x,t)=A_{\varepsilon}(x, t-{\rm{min}}(t,\tau_S))-A_{\varepsilon}(x,t)$ with a standard integral $A_{\varepsilon}(x,t)$ given in Appendix A.
Equation (\ref{sol}) holds for times $t$ between
$\max_{k,m}\{T_k^m\}$ (i.e., the last firing event from the set of all spine firing times), and  $T_{j}^\ell$
the time of the new firing event at the $\ell$th spine.


The firing times for the construction of solution (\ref{sol}) may be found from the set of threshold conditions $U_n(t)=h$, with $U_n(t)$ obtained by integrating equation (\ref{sp}).
In particular, to find a new firing time $T_{j}^\ell>\max\{T_k^m\}$ corresponding to the spine at location $x_{j}$ we have to solve the set of threshold conditions for the functions
\begin{equation}
U_n(t) = \frac{D r_a}{\widehat{C}r^2}\sum_{k,m} \widehat H(x_n-x_k,t-T_k^m)-h\sum_m  \e^{-\varepsilon_0 (t-T_n^m)},
\label{sol sp}
\end{equation}
where
\begin{equation}
\widehat H(x,t) = \int_{0}^t {\rm e}^{\varepsilon_0(s-t)} H(x,s){\rm d}s ,
\label{Hhat}
\end{equation}
and $\varepsilon_0=(1/\widehat{r}+1/r)/\widehat{C}$.
For $\varepsilon>\varepsilon_0$ this integral can be found in closed form as
$
\widehat H(x,t) = (A_{\varepsilon}(x,0)({\rm e}^{-\varepsilon_0(t-{\rm min}(t,\tau_S))}-{\rm e}^{-\varepsilon_0 t})+\widehat A(x,t-{\rm min}(t,\tau_S))-\widehat A(x,t))/\varepsilon_0
$
with 
\begin{equation}
\widehat A(x,t)={\rm e}^{-\varepsilon_0 t}\left[A_{\varepsilon-\varepsilon_0}(x,0)-A_{\varepsilon}(x,0)-A_{\varepsilon-\varepsilon_0}(x,t)\right]
+A_{\varepsilon}(x,t).
\label{eq:ahat}
\end{equation}
Alternatively the integral in (\ref{Hhat}) can be readily evaluated
numerically for the explicitly given function $H(x,t)$.  

By solving the set of threshold conditions with $U_n(t)$ defined by
(\ref{sol sp}) we obtain a vector of times showing
when each spine is able to reach the threshold
$h$. The smallest time from this vector, $T_{j}^\ell$, that satisfies
the refractory restriction $T_j^\ell-T_j^{\ell-1}>\tau_R$ defines a
new spiking event at  location $x_j$. 
The firing times of an individual spine have to be
separated by at least $\tau_S$, so the refractory time is
restricted by $\tau_R\geq \tau_S$.
As a result of finding the newly fired spine extra terms have to be added into both sums in
(\ref{sol sp}). The same routine is then repeated to obtain subsequent firing events.


\subsection{Numerical implementation of the SDS model}
\label{implementation}


The analytical integration of the equations of motion to get 
the explicit equations (\ref{sol}) and (\ref{sol sp}) for the partial SDS model
obviates the need for the numerical solution of a partial differential
equation. 
In the computations below we have taken
$\varepsilon>\varepsilon_0$ so that the explicit solution (\ref{Hhat})
and (\ref{eq:ahat}) can be used.
Up to determining the firing times this gives us an analytic solution.
The numerical scheme for the explicitly defined
solution was implemented in {\sc matlab}
and the firing times (defined from the threshold conditions $U_n(t)=h$) were
determined numerically using the root-finding routine {\tt fzero}.  
For each time interval the number of equations used for 
finding the earliest threshold crossing event
can be reduced by excluding the spines that are still in the refractory state.
Once the latest firing event has been determined then $V(x,t)$ can be evaluated (in terms of a sum over all spines and all previous firing events).  Since contributions from firing events in the distant past are exponentially decaying it would be safe to truncate this sum over events, and further improve numerical efficiency, though we do not do so in this paper.


In the following section we compare and validate the results from the (quasi)
analytic solution of (\ref{Hhat}) and (\ref{eq:ahat}) detailed above.
We show convergence by solving the explicit equation
(\ref{firstorder}) for $V(x,t)$ rather than (\ref{sol}) which requires
the solution of the differential equation (\ref{sp}) for each spine.
We also compare to the SDS model given as a system of differential equations
implemented in {\sc neuron} \cite{Hines}. We solve the partial model with
flux $I_{\rm sp}=\widehat V/r$ from the spine-head into the cable
as well as the full SDS model (\ref{cable}) and
(\ref{sp}). We also compare to the BR model with discrete
spines, again implemented in {\sc neuron}.  All code is available upon request to the authors.

\section{Validation of models and saltatory wave propagation}
\label{saltatory}

We validate the different models and their numerical solution by
comparing the propagation of a solitary saltatory wave. 
We start by noting that for the (quasi) analytic case a solitary wave 
can be determined from the solution of a nonlinear algebraic equation.

Assume that the spines are regularly distributed along the cable with spacing $d$, i.e., $x_n=nd$.
The speed of a solitary wave propagating in such a system may be determined in a self-consistent manner (see \opencite{Coombes03}) by the implicit equation
\be
h=\frac{\Lambda}{\widehat C r}\sum_{n=1}^{\infty}\widehat H(nd,n\Delta) ,
\label{continuumSDS}
\ee  
for the partial SDS model.
Here $\Delta$ measures the time between successive threshold crossings at adjacent spine-heads
and thus, the speed is $v =d/\Delta$. In Figure \ref{speed_xpp} we
plot this speed (solid curve) as a function of distance between the
spines. If the spines are separated beyond some critical value (to the
right of  limit point LP) the wave fails to propagate.  
From the analysis  in \cite{Coombes01a} it is further possible to establish that
the faster of the two branches is stable.
On the same figure we plot the results obtained from 
direct numerical simulations.
In all simulations the wave is initiated from an activated single spine at one end of the cable (and free boundary conditions are assumed).
Crosses indicate the 
speed obtained by using the explicit solution of the SDS model
given by equations (\ref{sol}) and (\ref{sol sp}).
Stars demonstrate the speed found by solving the system of
differential equations. There is  excellent agreement
between both approaches as well as with the speed found from equation
(\ref{continuumSDS}). Circles  in Figure \ref{speed_xpp} denote the
results of the full SDS model. The spine current passed into the cable
in the full model is relatively less than in the partial model 
when other parameters are fixed.
This explains why the wave in the last model propagates faster with
the LP for propagation failure shifted to the right. 
If the firing events in the model are defined by using the full form of the solution for $V(x,t)$ given by equation (\ref{fullsol}) the results will agree with the numerical solution of the full SDS model.
To demonstrate this convergence of the explicit approach, we evaluate
the speed of wave propagation 
using the explicit equation (\ref{firstorder}) for $V(x,t)$ rather than
(\ref{sol}).
Equation (\ref{firstorder}) includes a second order term in $\Lambda$ and this term has to be evaluated numerically.
In this case the functions $U_n(t)$ used for the determination of firing events are found by solving the differential equation (\ref{sp}) for each spine.
The obtained speed of the wave is indicated by triangles in Figure \ref{speed_xpp}. The inclusion of the second order term in the expansion for the full solution yields a slower wave speed (than with just the first order term). 
As more terms are included in the expansion the wave speed gets closer to the solution of the full model.
The inner plot in this figure is a magnified view (for small lattice spacing), 
and demonstrates that the speed of the wave attains a maximum for small $d$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{speed_SDS.eps}
\caption{A plot of solitary wave speed $v$ (in units
  $\lambda/\tau$) as a function of the distance between the spines $d$
  (in units $\lambda$). Solid line: the solution of equation
  (\ref{continuumSDS}). Crosses: the solution given by the explicit
  equations (\ref{sol}) and (\ref{sol sp}). Stars and circles: the
  solutions defined by solving (\ref{cable}) and (\ref{sp}) for the
  partial and full SDS model respectively. Triangles: the solution
  defined by solving (\ref{sp}) and (\ref{firstorder}). Parameters:
  $D=1$, $\tau_S=\varepsilon=1$, $\eta_0=1$, $\varepsilon_0 = 0.8$,
  $\widehat C=2.5$, $r=1$, $\tau_R=10$, $\widetilde{h}\equiv h/D
  r_a=0.05$. Inner plot is a magnified view of a part of the speed
  curve showing the LP for $v$.} 
\label{speed_xpp}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In Figure \ref{profile time} we plot the voltage at the location of the 10th spine along the cable when $d=0.4$. 
Here it is seen that the voltage of the fast wave in the partial SDS model has a higher amplitude compared to that of
the full SDS model.  However, even working to only second order in $\Lambda$ yields an improved solution, with increasing agreement as we include more terms in the series expansion for $V(x,t)$.

\subsection{Validation against the Baer-Rinzel model}

The SDS model defined in terms of the IF process is a
reduced version of the biophysically realistic model of  
\citeauthor{Baer91} where the spine-head dynamics are modeled by
HH kinetics. To compare these two models, we consider 
the BR model with a discrete distribution of spines
rather than with the spines uniformly distributed along the cable
(\opencite{Baer91}, \opencite{Lord02}).  In this case the model is defined by
equation (\ref{cable}) and
\be
\widehat C \PD{\widehat V}{t}=-I(\widehat V,m,n,h)- \frac{\widehat V-V}{r}.
\ee 
Here the function $I$ is the standard HH current, listed in Appendix B for completeness.
The membrane current arises mainly through
the conduction of sodium and potassium ions through voltage dependent
channels in the membrane.   A numerical simulation of the BR model with a
regular distribution of spines is presented in Figure
\ref{BR} which shows a space-time density plot of $V(x,t)$.
This plot nicely illustrates an example of a saltatory traveling wave
solution.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{BR.eps}
\caption{An example of saltatory traveling wave in the BR model with a discrete distribution of spines and spine spacing $d=5$. Other parameters are given in Appendix B.} 
\label{BR}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The speed of a saltatory wave changes when we vary the distance between the spines, and is quantified in Figure \ref{speed_BR}.   Also in this figure we show a plot of the wave speed for the partial SDS model (solid curve), with parameters chosen so as to give a corresponding fit to the more biophysically detailed BR model.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{SDS_BR2.eps}
\caption{The speed of wave propagation $v$ as a function of the distance between the spines $d$ in the discrete BR model (dotted curve) and in the partial SDS model (solid curve). Parameters of the SDS model:
$D=0.0003$ m/s, $\eta_0=100$ mV, $\tau_S=2$ ms, $r=0.01$ M$\Omega$, $\widehat C=0.001~\mu$F, $\varepsilon=3$, $\varepsilon_0=0.7$, $\tilde{h}\equiv h/r_a=0.00095$. Parameters of the BR model are given in Appendix B.} 
\label{speed_BR}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

From now on we consider the SDS model to be the (quasi)
analytic solution of the partial SDS model obtained by solving equations
(\ref{sol}) and (\ref{sol sp}) with a numerical approximation to the
firing time. 

\subsection{Saltatory Waves and spine distributions}

Waves generated in the SDS model travel in a saltatory manner 
because the discreteness of the spine distribution breaks the translation symmetry of the underlying cable.
An example of a saltatory wave (constructed from (\ref{Hhat}) and (\ref{eq:ahat})) is shown in
Figure \ref{profile space}.  
In Figure \ref{pulse} we show a space-time plot of the cable voltage
with 30 regularly distributed spines, which nicely illustrates (as does Figure \ref{BR}) that although
the wave has a saltatory nature it  travels with well-defined speed. 
For sufficiently small values of the refractory period it is also possible to generate periodic traveling waves.
An example is shown in Figure \ref{periodic}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{d04full_part_mid.eps}
\caption{Voltage $V$ (in units of $\eta_0$) at the location of the 10th spine along the cable  with spine spacing $d=0.4$. Other parameters as in Figure \ref{speed_xpp}. 
Left (right) solid curve: the solution of the partial (full) SDS model. Dashed curve: the solution of the full SDS model including terms up to order $\Lambda^2$.
}
\label{profile time}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{profile_d.eps}
\caption{An example of saltatory solitary wave in the partial SDS model traveling  along the cable with spine spacing $d=0.6$ and $D=0.25$.
Other parameters as in Figure \ref{speed_xpp}.
Snapshots are shown for times $t_0=7\Delta+\tau_S$, $t_1=8\Delta+0.1\tau_S$, $t_2=8\Delta+\tau_S$ and $t_3=9\Delta+0.1\tau_S$, where $\Delta=1.1217$.}
\label{profile space}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{pulse_last.eps}
\caption{An example of saltatory traveling pulse in the partial SDS model initiated from a
  single active spine and moving out through the cable with 30
  regularly spaced spines when $d=0.6$. Other parameters as in Figure
  \ref{speed_xpp} except $\tau_R=6$.} 
\label{pulse}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.5in]{periodic.eps}%{irreg_pulse_4e.eps}
\caption{An example of periodic traveling wave 
for the model of Figure \ref{pulse} with $\tau_R=5$.
Plots on the right illustrate the voltage potential as a function of time $t$ at the location of the 10th spine along the cable (top plot) and as a function of space $x$ at time $t=25$ (bottom plot).}
\label{periodic}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In real dendritic tissue  the distribution of spines 
is likely to be irregular. 
In the SDS framework we might simply consider that the distribution function $\rho(x)$ has some disorder.
An example of the effect that such disorder can have is shown in Figure \ref{irreg}, where we choose spine positions from a uniform distribution.
In comparison to Figure \ref{pulse}, not only does this cause irregular (and repeated) wave propagation, but it may lead to back-propagating waves.  As expected, too much spatial disorder in the spine distribution  can lead to propagation failure, when any one pair of adjacent spines becomes too far separated to cause firing of one by the other. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=3.5in]{back_prop1.eps}%{irreg_pulse_4e.eps}
\caption{An example of irregular wave propagation in the partial SDS model for the parameter set corresponding to Figure \ref{pulse} with 33 spine positions chosen from a uniform distribution.}
\label{irreg}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Filtering properties}
\label{filtering}

The spatio-temporal patterns of activity in neurons and indeed networks of neurons is intimately linked 
to the processing of sensory information.  A fundamental issue in understanding  the relationship between structure and function is how neurons transform and represent information using such patterns.  
In this regard an understanding of how neural systems respond to sensory stimuli is an important issue.  Indeed, although the selective behavioral responses to stimuli in the sensory systems of many organisms is believed to rely on the presence of filters for the neural representations of the temporal structure of sensory signals, the mechanisms underlying this temporal selectivity are not completely understood.  Undoubtedly these mechanisms are subserved by the biophysics of neurons, 
including the passive electrical properties of a neuron, and the types and distribution of channels and conductances they support.  However, only relatively recently has a role for dendritic spines begun to be seriously entertained.
For an increased understanding of the role of spines in temporal filtering, intracellular recordings are necessary. 
Progress in recording from neurons as small as 10 $\mu$m has been made using `patch-type' pipette techniques \cite{Rose96},  and these have been used to investigate the role of spines in filtering temporal input patterns 
to midbrain neurons of the weakly electric fish {\it Eigenmannia} \cite{Fortune97,Rose99}.
In particular, this work has shown that neurons with a  broad dendritic arbor and relatively spiny dendrites demonstrate low-pass temporal filtering properties in response to input stimuli.
Here, we suggest the use of the SDS model for re-examining the findings of Rose and Fortune from a theoretical perspective, with a particular emphasis on clarifying the effect of both passive and active properties of the dendrite on filtering.
The SDS model incorporates passive membrane properties within the cable and mimics active properties by the IF process in the spine-head.  Importantly, we have the ability to vary spine distributions and see the effect on filtering.  A limitation of the model is that no explicit consideration is given to the dynamics in the soma.
However, output signal recordings at the soma in experiments are believed to be primarily influenced by the spread of activity  in the dendrites, suggesting that an analysis of filtering within the SDS model will still be relevant.

We assume that an electrical signal in the form of a periodic pulse train is applied to the dendritic cable at position $x_0$. 
This is modeled by adding the following term to the right-hand side of equation (\ref{cable1})
\be
I_{\rm{stim}}=\sum_p\delta(x-x_0)\delta(t-pT) ,
\label{Istim}
\ee
where $T$ is the period of the stimulus.
By solving equation (\ref{cable1}) with this applied current, we obtain an extra term that enters  the right hand side of equation (\ref{sol}) for the membrane potential as follows
\be
\sum_{p=0}^P G(x-x_0,t-pT),
\ee 
where $P=\max \{p~|~t-pT\geq 0\}$.
Therefore, equation (\ref{sol sp}) for the spines has to be updated accordingly by adding an extra term
\be
\frac{1}{r \widehat C} 
\sum_{p=0}^P \widehat{G}(x_n-x_0,t-pT),
\ee
where $\widehat G(x,t)$ is found as 
\be
\widehat G(x,t)=\int_{0}^t G(x,s)\e^{\varepsilon_0 (s-t)} {\rm d} s = 
\e^{-\varepsilon_0 t}(A_{\varepsilon-\varepsilon_0}(x,0)-A_{\varepsilon-\varepsilon_0}(x,t)). 
\ee


In Figure \ref{d04T_all} we plot the results of simulations of the SDS model driven by the applied forcing signal (\ref{Istim}). In all simulations
the distribution of spines is chosen to be regular and the current is applied to the left end of the cable at a distance of $0.5$ (in units of $\lambda$) from the first spine. In each of these three examples we stimulated the system with a different fixed frequency. Differing applied stimuli lead to different patterns of waves propagating through the system. The signal transformation during its propagation is best established by comparing the applied stimulus and the system response at the opposite end of the cable.
Plots on the right in Figure \ref{d04T_all} show the membrane voltage at the location of the $55$th spine along the cable.
If the period of stimulation is large, the signal measured at the end of the cable has a single inter-spike interval (ISI) consistent with the period of stimulus (see Figure \ref{d04T_all}A). If we decrease the period of stimulation making it closer to the refractory time $\tau_R$ the output response at the end of the cable shows the presence of two ISIs (in Figure \ref{d04T_all}B).
When the period of applied current is chosen to be small
(Figure \ref{d04T_all}C),
the refractoriness of the system prevents the output from  reaching the frequency of the stimulus. 
Instead, the system initiates  repetitive waves consistent with the refractory time, with a single ISI. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.5in]{d04T150_60_10.eps}
\caption{Examples of wave propagation in the SDS model driven by periodic stimulation.
The current stimulus is applied to the left end of the cable 
embedded with 60 regularly spaced spines.
The system parameters as in Figure \ref{speed_xpp} with $d=0.4$, except $\tau_R=7$.
The period of stimulation is $T=15$ (A), $T=6$ (B) and $T=1$ (C).
Plots on the right are examples of voltage profile at the location of the 55{\it {th}} spine along the cable.}
\label{d04T_all}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The observed correlation between the periods of the applied signal and the ISIs of the system responses  is quantified in Figure \ref{freq}.
The $x$-axis of the plots in this figure contains the frequencies of the applied stimulus used in the simulations.  Along the $y$-axis we plot the frequencies at the end of the cable measured as the number of spikes emitted per unit time. Figure \ref{freq}A shows the data from simulations for a spine spacing with $d=0.4$. Vertical dashed lines in this figure indicate the frequencies of the applied signal used in simulations in Figure \ref{d04T_all}. The two values for the spike frequencies
that are defined for the particular range of input frequencies illustrate the presence of two distinct ISIs, as  for example in Figure \ref{d04T_all}B. 
The response of the system coincides with the input signal for low frequencies of stimulation, whereas for high frequencies  it is limited by the system refractoriness. These dynamics indicate that the SDS model exhibits  low-pass temporal filtering properties. 
The doubly periodic ISIs occur at stimulus periods close to the refractory time-scale. The smaller the distance between the spines, the smaller the range of input frequencies that lead to doubly periodic ISIs. However, if the spines are further apart, the regime of stimulus frequencies where the system responds irregularly with doubly periodic ISIs is increased.
An example is shown in Figure \ref{freq}B, where the spine spacing $d$ equals 0.8. 
Thus, the low-pass filtering properties are reduced with decreasing spine density. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.5in]{Freq.eps}
\caption{Frequency of the applied stimulation vs 
frequency of spikes in the output response at the end of the cable. The system parameters as in Figure \ref{d04T_all} with $d=0.4$ in A and $d=0.8$ in B.
Vertical dashed lines in A indicate the frequencies
of stimulation that were chosen for simulations in Figure \ref{d04T_all}.}
\label{freq}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In Figure \ref{amp} we show the relative amplitude of voltage responses to the frequency of input current. This figure is included to make a link to the experimental observations that spiny neurons demonstrate a decline in the amplitude of voltage responses with increasing stimulation frequency.
The relative amplitude measured in dB is evaluated as the difference between the maximum amplitude and the mean of the signal response.
Figure \ref{amp}A is plotted for two different spine spacings $d$. A decrease in the spine density (achieved by increasing $d$) reduces the decline in the voltage amplitude. 
The presence of doubly periodic ISIs in the voltage response causes the  irregularities seen in both curves 
at periods around the refractory time-scale.
In Figure \ref{amp}B we show that 
a change in the refractory time in the system results in a shift of the original curve but does not affect the overall trend.

To summarize, the filtering properties of the SDS model are entirely consistent with the experimental observations of Rose and Fortune \cite{Fortune97,Rose99}, namely that spiny dendrites with active membrane properties demonstrate low-pass temporal filtering properties.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.50in]{amp.eps}
\caption{Frequency of applied stimulation vs relative amplitude of the voltage potential (in dB) 
with the system parameters as in Figure \ref{d04T_all}.
A: $d=0.4$ (squares), $d=0.8$ (triangles), $\tau_R=7$.
B: $\tau_R=7$ (squares), $\tau_R=9$ (triangles), $d=0.4$.}
\label{amp}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Noise induced phenomenon}
\label{noise}


The spontaneous behavior of neurons \textit{in vivo} is believed to be driven by  voltage fluctuations arising from system noise. 
Such noise can be characterized as either intrinsic or extrinsic.
A major source of intrinsic noise arises from the stochastic gating of ion channels in the cell membrane.
Importantly, if the membrane potential is close to threshold, channel noise can be critical for the generation of action potentials. In the SDS model the description of noise due to stochastic channels is best introduced  
in the excitable spine-heads. However, since
the kinetics of ion channels are not explicitly included in the model, the noise is simply considered to lead to the stochastic generation of action potentials.
Extrinsic noise sources in the neuron typically arise from synaptic inputs.  Interestingly, there is an observed difference  between the spontaneous activity of neurons {\it in vivo} and the activity during intracellular stimulation {\it in vitro}. 
Introducing a voltage fluctuation in the cable membrane as well as in the spines is a natural way to mimic the presence of input noise in the SDS model of an isolated neuron.

To study the effects of both types of noise on the properties of wave propagation in the SDS framework we assume that the system is driven by  additive
white noise in time that is either smooth (correlated) or white
(uncorrelated) in space. 
The SDS model that includes noise both in the cable and the spine-head
is given by 
\ba
\label{cable noise}
 {\rm d} V &=& \left[ D \PD{^2 V}{x^2}-\frac{V}{\tau}+D
 r_a\rho(x)\frac{\widehat{V}-V}{r}\right] {\rm d} t +\mu_V {\rm d} W_V(t,x),\\ 
{\rm d} U_n &=& \left[ \frac{V_n}{\widehat{C} r}-\varepsilon_0
U_n-h\sum_m\delta(t-T_n^m)\right] {\rm d}t +\mu_U {\rm d} W_U(t,x). 
\label{sp noise}
\ea
This is in fact a system of integral equations 
where we interpret noise in the It\^o sense.
We assume that the Wiener processes $W_{V,U}$ are of the form 
$$W(t,x)\,=\,\sum_{n\in\mathbb{Z}}b_n\phi _n \beta
_n(t)~, $$
where the $\beta_n$ are a mutually independent ordinary set of Brownian motions and the 
$\phi_n$ are appropriate basis functions (with either periodic or
Neumann boundary conditions). The coefficients 
$b_n$ determine the spatial correlation. 
For spatially correlated noise we take 
$$
\sum_{n\in\mathbb{Z}}e^{2\alpha |n|}|b_n|^2\,<\,\infty~.
$$
Here, the parameter $\alpha$ determines the correlation length scale and the 
spatial smoothness. 
Examples of such noise are considered for example in
\cite{DaPrtoZbczyk,Shardlow05,GrcaOjlvo+Sncho,Lord03}.
The parameters $\mu_V$ and $\mu_U$ describe the strength of the noise in the
cable and in the spine-head respectively.

The integration of equation (\ref{cable noise}) for
the membrane potential  leads to a solution for $V(x,t)$ that may be split into two terms.
The first term, representing the deterministic component, may be
computed using the exact solution given by equation (\ref{sol}),
whereas the second term incorporating the noise has to be evaluated
numerically.  To evaluate this term we enforce a spatial correlation in
Fourier space and use a Brownian increment of mean zero and finite variance in time, as discussed for example in \cite{kloeden92}. 
The noise for an infinite cable was approximated by applying periodic boundary
conditions of twice the numerical cable length
and taking the central
portion away from the boundary.
The integration of equation (\ref{sp noise}) for $U_n(t)$, similarly to
the cable, generates two terms.  
In general, the first term has to be evaluated numerically and we use
an explicit Euler-Maruyama method of approximation.  The second term,
that defines the noise in the spine-head, was generated in the same
manner as for the cable and then sampled at the spine position.
However,  when $\mu_V=0$, the model reduces to the case where only
the spines are forced by noise and, thus, the solution for
$V(x,t)$ is explicitly defined. In this limit, the numerical
evaluation of the  first term in the solution of equation (\ref{sp
  noise}) can be avoided. Instead, this term is found using the
explicit equation (\ref{sol sp}). 


Here we demonstrate the effect noise can have on the properties
of wave propagation in the SDS model. 
We begin by exploring the SDS model with noise only in the
spine-head and take $\mu_V=0$. 
In all of the figures that follow the noise-path is the same and
$\mu_U$ is the only parameter that varies.
In all simulations the spines are regularly distributed along the
cable.  
When $d=0.6$ (left of LP in \ref{speed_xpp}) there is propagation in the absence of noise 
as illustrated in Figure \ref{sp wh d0_6}.
In general, by increasing the level of noise in the system
the patterns of cell response change from isolated or repetitive wave
propagation to the case of almost simultaneous firing of all
spines (i.e., a coherent behavior). However, as shown in Figures \ref{sp wh d0_6}A and B, at 
low noise levels the repetitive propagation of waves can be
suppressed by a small noise increase. 
For high noise (Figure \ref{sp wh d0_6}D) the majority of spines
can generate an action potential immediately and this leads to
the activation of residual spines (i.e., those that did
not fire) shortly thereafter. 
When $d=1$ (right of LP in \ref{speed_xpp}) the system does not support
traveling waves in the absence of noise.
However, the examples in  Figure \ref{sp wh d1} illustrate that  it is possible
for noise to sustain spatio-temporal structures that could not
otherwise occur. In particular, for a low level of noise we observe
purely noise induced waves such as shown in Figure \ref{sp wh d1}A. 
Propagation  failure occurs if the level of noise is below some critical value. 
In regions of high noise it is possible to see coherent oscillations
in the system (Figure \ref{sp wh d1}B) where almost all spines
simultaneously generate action potentials. 
Choosing correlated noise in the spine-head instead of  white
noise does not yield a qualitative change in system
behavior. This is illustrated in Figure \ref{sp cor d1} where we plot
two examples of spatio-temporal structures generated by correlated
noise. The patterns of wave propagation are similar to the
examples shown in Figure \ref{sp wh d1}.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{sp_wh_d0_6.eps}
\caption{Traveling waves in the SDS model in the presence
of uncorrelated white noise in the spine-heads for parameters $d=0.6$, $\mu_V=0$ and A: $\mu_U=0.007$, B: $\mu_U=0.009$, C: $\mu_U=0.015$, D: $\mu_U=0.028$. Other parameters as in Figure \ref{pulse} }
\label{sp wh d0_6}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{sp_wh_d1.eps}
\caption{Traveling waves in the SDS model in the presence
of uncorrelated white noise in the spine-heads for parameters $d=1$, $\mu_V=0$ and 
A: $\mu_U=0.03$, B: $\mu_U=0.055$.
Other parameters as in Figure \ref{pulse}.}
\label{sp wh d1}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{sp_cor_d1.eps}
\caption{Traveling waves in the SDS model in the presence
of correlated noise in the spine-heads for the parameters
$\alpha=1$, $d=1$, $\mu_V=0$ and A: $\mu_U=0.17$, B: $\mu_U=0.4$.
Other parameters as in Figure \ref{pulse}.}
\label{sp cor d1}
\end{center}
\end{figure}

The noise term in equation (\ref{sp noise}) for the spine-heads dynamics
leads to a stochastic process for spike generation.
Another way to incorporate this stochasticity 
into the firing events is to introduce a source of noise at the threshold level. This can be modeled under the replacement $h \rightarrow h+\xi$ where $\xi$ is an additive noise with distribution $\rho(\xi)$. This approach can be simply implemented by generating a random vector ${\bf \xi}$ of length $n$ from a distribution of mean $h$ and standard deviation $\beta$ at each small time step. Then the firing events are defined by checking the threshold conditions $U_n>\xi_n$ where $\xi_n$ is the $n$th element of vector $\xi$ and $U_n$ is given by equation (\ref{sol sp}). 
The level of noise in the system is controlled by the parameter $\beta$. The numerical simulations of the SDS model in the presence of
threshold noise demonstrate  consistency with the results of the model driven by the stochastic forcing at the spine-heads. 
It is also possible to model the effect of threshold noise using a 
probabilistic rule for spike generation.  To see this we consider drawing the threshold noise from a distribution $\sigma$, so that
the probability of a firing event can be written
\be
P(U_n>h)=\int\sigma(\xi)\Theta(U_n-h-\xi){\rm d} \xi =f(U_n-h),
\ee
where $f(\xi)=\int^\xi \sigma(x) {\rm d} x$.  Considering a bell-shaped noise distribution for $\sigma$, then the function $f$ will have a sigmoidal shape.
Here we use the following form for the function $f$
\be
f(U)=(1+{\rm e}^{-\gamma h})/(1+{\rm e}^{-\gamma U})-{\rm e}^{-\gamma h},
\label{prob}
\ee
so that the probability of a firing event is zero and one respectively for $U=0$ and $U\rightarrow \infty$.
The parameter $\gamma$ in (\ref{prob}) controls the level of noise in the system
by determining the width of the bell-shaped distribution $\sigma$ (such that decreasing $\gamma$ corresponds to increasing noise). 
Simulation results of the SDS model 
with this probabilistic update rule are shown in Figure \ref{sigmoid}. Plots on the left are for the system with spine spacing $d=0.6$, whereas plots on the right are obtained when $d=1$. For low noise the system is able to support repetitive wave propagation as illustrated in Figures \ref{sigmoid}A and D. When the noise level is below some critical value, propagation failure occurs. This is similar to the behavior observed in our earlier examples of the stochastic SDS model. Here, however, an 
increase in noise causes more irregularity in the patterns of wave propagation (see Figures \ref{sigmoid}B and E). Moreover, high noise limits the spread of activity in the system by inhibiting  wave propagation (see Figures \ref{sigmoid}C and F) or can even lead to propagation failure. 
The probabilistic rule used for the determination of the spikes in these last examples ensures that the probability of an individual spine being fired is low if the membrane voltage in the cable at the location of this spine is low. 
This explains why  wave propagation can terminate for high noise levels. 
In previous examples, where the noise in the spines is modeled by the stochastic differential equation (\ref{sp noise}), the high noise in the system, contrary to this last example, stimulates activity along the whole cable.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.5in]{sigmoidal_noise1.eps}
\caption{Traveling waves in the SDS model in the presence of threshold noise incorporated by the probabilistic rule for the spine spacing $d=0.6$ (A, B and C) and $d=1$ (D, E and F). The noise levels are $\gamma=60$ (low noise) (A and D), $\gamma=5$ (medium noise) (B and E), $\gamma=0.8$ (C) and $\gamma=1$ (high noise) (F). 
Other parameters as in Figure \ref{pulse}.}
\label{sigmoid}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now consider the stochastic SDS model driven by noise in the cable and thus, take $\mu_U=0$. In Figure \ref{cab wh d06} we plot results of simulations
for the system in the presence of uncorrelated white noise.
These examples demonstrate that in the parameter
region where waves propagate in the absence of noise ($d=0.6$),
an increase in noise level can lead to repetitive wave initiation. Plots on the right show profiles of action potentials in the membrane at the locations of the $20$th and $12$th spine along the cable (Figures \ref{cab wh d06}A and B).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{cab_wh_d0_6.eps}
\caption{Traveling waves in the SDS model in the presence
of uncorrelated white noise in the cable for the parameters
$d=0.6$, $\mu_U=0$ and A: $\mu_V=0.08$, B: $\mu_V=0.114$.
Other parameters as in Figure \ref{pulse}. Plots on the right show the profiles of membrane voltage at the locations of $20$th (A) and $12$th spine (B) along the cable.}
\label{cab wh d06}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the presence of correlated rather than white noise, the system demonstrates a similar kind of behavior.
This is shown in Figure \ref{cab cor d06} where increasing noise levels generate a single, double and repetitive wave respectively.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{cab_cor_d0_6.eps}
\caption{Traveling waves in the SDS model in the presence
of correlated noise in the cable for the parameters
$\alpha=1$, $d=0.6$, $\mu_U=0$ and A: $\mu_V=0.4$, B: $\mu_V=0.8$ and C: $\mu_V=0.81$.
Other parameters as in Figure \ref{pulse}. Plots on the right show the profiles of membrane voltage at the locations of $20$th spine along the cable.}
\label{cab cor d06}
\end{center}
\end{figure} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

If the distance between the spines is chosen to be large,  the noise in the cable, white or correlated, is not able to initiate any
wave propagation in the system. 
However, in the case of high white noise some limited activity might be possible.  An example is illustrated in Figure \ref{cab wh d1}. The plot on the right 
shows the firing times of spines that generated action potentials during the simulation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{cab_wh_d1.eps}
\caption{Traveling waves in the SDS model in the presence
of uncorrelated white noise in the cable for the parameters
$d=1$, $\mu_U=0$ and $\mu_V=0.35$. Other parameters as in Figure \ref{pulse}. Plot on the right shows the spines that generate the action potential.}
\label{cab wh d1}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


In summary, from our numerical simulation studies we conclude that 
noise in the spine-heads has a strong effect
on the properties of wave propagation. In the parameter regime where
the deterministic model fails to support waves,
noise in the spines may aid in the initiation and propagation of waves.
Moreover, for high noise in the spines one observes coherent oscillations. 
At the same time, the
stochastic SDS model shows robustness of wave propagation  to
low noise in the spine-heads. On the other hand, 
noise in the cable does not have a strong effect on 
patterns of activity.

We now consider the presence of noise both in the cable and in the
spine-heads (see  Figure \ref{noise cab sp}).
Noise in the cable was chosen as for the model in Figure 
\ref{cab cor d06}B. For Figure \ref{noise cab sp}A (Figure 
\ref{noise cab sp}B)  noise in the spine-heads was generated as in
Figure \ref{sp wh d0_6}B (Figure \ref{sp wh d1}A). Simulations
demonstrate that in the parameter regime where the reduction of noise to zero would still support
wave propagation the
cable noise stimulates the repetitive initiation of waves (see Figure
\ref{noise cab sp}A).  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.8in]{noise_cab_sp.eps}
\caption{Traveling waves in the SDS model in the presence
of uncorrelated white noise in the spines-heads and correlated noise in the cable for the parameters
A: $d=0.6$, $\alpha=1$, $\mu_V=0.8$, $\mu_U=0.009$ and
B: $d=1$, $\alpha=1$, $\mu_V=0.8$, $\mu_U=0.03$.
Other parameters as in Figure \ref{pulse}.}
\label{noise cab sp}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Finally, we investigate the robustness of the filtering properties observed in the deterministic SDS model to noise sources. The SDS model in the presence of a small amount of correlated noise both in the cable and the spine-heads was driven by injected current in the form of a periodic pulse train as in section \ref{filtering}. 
In Figure \ref{Ts} we demonstrate examples of membrane voltage profiles along the cable for three different periods of stimulation. These plots should be compared with the results of simulations shown  
in Figure \ref{d04T_all}.
The presence of noise in the system for the low and high periods of stimulation (Figures \ref{Ts}A and C) retains the same ISIs as in the deterministic model. When the period of stimulation is close to the refractory time the well-defined pattern of doubly periodic ISIs seen in Figure \ref{d04T_all}B can be destroyed by noise. However, in general the low-pass temporal filtering properties observed in the deterministic model are robust to noise.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=4.5in]{Ts.eps}
\caption{Examples of voltage profiles at the location of the $55$th spine along the cable in the stochastic SDS model driven by the applied current stimulation of period $T$. The parameters are $\alpha=0.4$, $\mu_V=0.49$, $\mu_U=0.043$ and $T=15$ (A), $T=6$ (B) and $T=1$ (C). Other system parameters as in Figure \ref{d04T_all}.}
\label{Ts}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




\section{Discussion}
\label{discussion}

It is now  just over one hundred years since the discovery of spines, yet there are still mysteries about the contribution they make to single neuron dynamics.  Building on the insight into dendritic function gained since the pioneering  theoretical work of Rall in the 1950's (surveyed in \cite{Segev95}), we propose that the SDS model is ideally suited to probing how active spines shape and modulate a spatially structured input.  In this paper we have focused on a rather simple (spatially localized) input, though we emphasize here that more spatially structured input can also be treated.  Moreover, there is no barrier to working with a truly branched dendritic tree model as the analytical (``sum-over-paths'') techniques necessary to do this have already been developed by Bressloff et. al. (\citeyear{Bressloff96}).
In future work we propose to use the SDS model to understand whether inputs to a neuron sum linearly, or depend on spatial relationships in the input, such as clustering.  The currently available experimental evidence is conflicting \cite{Cash98,Polsky04}.  
Interestingly, recent advances in imaging technology (using fluorescent dyes in combination with confocal or two-photon laser scanning microscopy) have permitted time-lapse observation of spine morphology in living neurons \cite{Fischer98,Dunaevsky99,Bonhoeffer02}, showing that spines are also constantly moving.  In fact over time-scales of seconds, spines continuously undergo small changes in shape, thought to be powered by dynamic actin filaments. On timescales of minutes to hours spines can change their shape dramatically or even appear or disappear. Moreover, new spines can be generated in response to synaptic stimulation that also results in strengthening of synapses. These morphological changes could underlie some of the changes in synaptic strength induced by neural activity.  
Working within the SDS framework we would also be in a position to examine how successive synaptic input might influence spine motility and electro-chemical properties.  Some work in this direction has already been done by Verzi et. al. (\citeyear{Verzi05}), for the continuum Baer and Rinzel model using phenomenological models of spine density dynamics.  Within the SDS framework we envisage complementary work, building on existing modeling studies (of spines as calcium compartments) such as that of Holcman and Schuss (\citeyear{Holcman05}), exploring how calcium accumulates in discrete spines, and how this makes its way to the cell nucleus and triggers genetic mechanisms that ultimately lead to the mechanical reconfiguration of the synapse.  Both the above issues are topics of current investigation and will be reported upon elsewhere.




\section*{Acknowledgments}
We would like to thank Steve Cox, Kre\v{s}imir Josi\'{c} and Costa Colbert for useful discussions about potential applications of the SDS model.  The work in this paper is supported through EPSRC Grant No. GR/S60914/01.
SC would also like to acknowledge ongoing support from the EPSRC through the award of an Advanced Research Fellowship, Grant No. GR/R76219.


\section*{Appendix A: $A_{\varepsilon}(x,t)$}

\begin{eqnarray}
\nonumber
A_{\varepsilon}(x,t)&=&\frac{\eta_0}{4}\sqrt{\frac{1}{\varepsilon D}}\left\{ {\rm exp}\left( -\left|x\right| \sqrt{\frac{\varepsilon}{D}} \right){\rm erfc}\left(-\frac{\left|x\right|}{\sqrt{4Dt}}+\sqrt{\varepsilon t}\right) \right. \\
&+&\left. {\rm exp}\left( \left|x\right| \sqrt{\frac{\varepsilon}{D}} \right){\rm erfc}\left(\frac{\left|x\right|}{\sqrt{4Dt}}+\sqrt{\varepsilon t}\right)  \right \}.
\nonumber
\end{eqnarray}
\section*{Appendix B:  Equations of the BR model}


\ba
\nonumber
I(\widehat V,m,n,h)&=&g_Kn^4(\widehat V-v_K)+g_{Na}hm^3(\widehat V-V_{Na})+g_L(\widehat V-V_L),\\
\tau_{X}(\widehat V) \FD{X}{t} &=& X_{\infty}(\widehat V)-X \nonumber,
\ea
for $X \in \{m,n,h\}$ where $V_L$, $V_K$ and $g_{Na}$ represent the constant membrane reversal potentials associated with the leakage, potassium and sodium channels respectively.
\ba
\nonumber
\tau_X(\widehat V)=\frac{1}{\alpha_X(\widehat V)+\beta_X(\widehat V)},
\ea
\ba
X_{\infty}(\widehat V)=\alpha_X(\widehat V)\tau_X(\widehat V) ,
\nonumber
\ea
where
\ba
\nonumber
\alpha_m(\widehat V)=\frac{0.1(\widehat V+40)}{1-{\rm exp}[-0.1(\widehat V+40)]},
\ea
\ba
\nonumber
\alpha_h(\widehat V)=0.07{\rm exp}[-0.05(\widehat V+65)],
\ea
\ba
\nonumber
\alpha_n(\widehat V)=\frac{0.01(\widehat V+55)}{1-{\rm exp}[-0.1(\widehat V+55)]},
\ea
\ba
\nonumber
\beta_m(\widehat V)=4.0{\rm exp}[-0.0556(\widehat V+65)],
\ea
\ba
\nonumber
\beta_h(\widehat V)=\frac{1}{1+{\rm exp}[-0.1(\widehat V+35)]},
\ea
\ba
\nonumber
\beta_n(\widehat V)=0.125{\rm exp}[-0.0125(\widehat V+65)].
\ea
The following parameter values were used: $R_a=36$ $\Omega \cdot$cm, $R_m=3333$ $\Omega \cdot$cm$^2$, $C_m=1$ $\mu$F/cm$^2$, $a=1$ $\mu$m, $r=0.05$ M$\Omega$, $\widehat C=0.001$ $\mu$F, $g_L=0.0003$ mS, $g_K=0.036$ mS, $g_{Na}=0.12$ mS, $V_L=-54.402$ mV, $V_K=-77$ mV and $V_{Na}=50$ mV.




\bibliographystyle{klunamed}
\bibliography{SDSfiltering}

\end{article}
\end{document}
