\graphicspath{ {./images/} }

\section{Biology background}
To understand the function and behavior of ion channels, several points need to
be discussed. The background
comes from a wide disciplinary range, including biology, chemistry,
thermodynamics, electricity, probability theory, and differential equations.
The most important of these points are presented in the following subsections.
For more information, the reader is referred to Hille's extensive book,
\emph{Ion Channels of Excitable Membranes} \cite{hille-2001}, and Keener and
Sneyd's work on the applied mathematics of cellular physiology
\cite{keener-2009}.


\subsection{Salts, ions, and the cellular solutions}
Before discussing the biology of cells and ion channels any further, a quick
discussion regarding ionic solutions and the
\emph{electrochemical gradient} is needed. 
As any good student should know, ions are atoms or molecules that have a
positive or negative charge due
to differing numbers of protons and electrons. Specifically, anions have more electrons than
protons, so they have a negative charge. Conversely, cations
are ``missing'' electrons, and therefore have a positive charge. 

Ions arise from a variety of natural processes, including the dissolution of salts
in water. For example, when table salt, or
more officially, sodium chloride
($NaCl$) dissolves in water, sodium cations with one positive charge ($Na^+$)
and chloride anions with one negative charge ($Cl^-$) become dissociated and
free to move around independently. Similarly, the salts potassium chloride
($KCl$) and calcium
chloride ($CaCl_2$) dissolve into solutions of free $K^+$, $Cl^-$, and $Ca^{2+}$
ions. 

Solutions of dissolved salts are of critical importance to a cell's function.
Many intracellular and extracellular environments contain (among other 
molecules) these four ions in very specific
concentrations. Ionic concentration values for three different cell types are provided 
in Table \ref{table:concentrations}.

Now let us discuss the concept of the electrochemical gradient. 
Consider the situation presented in Figure \ref{fig:nernst_potential} in which
a water-filled bath is separated into two halves by a \emph{semi-permeable} membrane.
This membrane contains small pores which allow $K^+$ ions, but nothing else, to 
pass through. Therefore, $K^+$ ions can move freely across the membrane from the left
to right side of the bath, and vise versa. A more concise description of the membrane 
is to say that it is perfectly permeable to $K^+$. 
As a convention, the left side of the bath represents the inside of a biological cell
(intracellular), while the right side is the extracellular environment.

%%%% Nernst potential  %%%%
\begin{figure}[h!]
\begin{center} 
\includegraphics[width=6in]
{nernst_potential.eps}
    \caption{A voltmeter measures the potential difference between the two sides
of a bath which is separated by a membrane perfectly permeable to $K^+$. The
left side contains a higher concentration of dissolved $KCl$ than the right. At
the first instant the solutions are added, no ions have had time to diffuse through the
membrane, so the charges are balanced resulting in a
voltage of zero. However, $K^+$ ions immediately begin diffusing through the
membrane, down their concentration gradient resulting in a buildup of positive
charge of the right side of the membrane. Eventually the electrical and
diffusion forces will balance and the system will reach equilibrium. The
equilibrium potential is equal to the Nernst potential of potassium, $E_K$.}
    \label{fig:nernst_potential}
\end{center}
\end{figure}

At time $t=0$, a concentrated solution of $KCl$ is introduced to the left side of the 
bath, while a weaker solution of $KCl$ is introduced to the right side. This configuration
is represented by the larger letters of $K^+$ and $Cl^-$ on the left side. 
A voltmeter measures the voltage between the left and right sides, and reads zero volts
at $t=0$, because the number of cations and anions are equal on both sides of the membrane.
Because the left side has a higher concentration of $K^+$, the ions experience a diffusion 
force which is represented by arrows crossing the membrane. $Cl^-$ ions, however, cannot 
cross the membrane as it is not permeable to them.

The instant after the solutions are added to the bath, $K^+$ will begin moving across
the membrane from left to the right sides, driven by the diffusion force. As this movement
of charge continues, the right side will become more positive. Eventually the buildup of positive
charge will deter $K^+$ ions from moving to the right side. At this point, the system
is at equilibrium. The chemical gradient pushing $K^+$ to the right is balanced by
the electrical force opposing the movement of charge. The voltage at which the
bath reaches equilibrium is called the Nernst potential, and is calculated for
a generic ion species $S$ like so:

\begin{equation}
 E_S = \frac{RT}{zF}ln \left(\frac{[S]_e}{[S]_i} \right),
 \label{eq:nernst}
\end{equation}

where $R$ is the gas constant, $T$ is the absolute temperature in
Kelvins, $z$ is the charge on ion $S$, $F$ is Faraday's constant, and
$[S]_e$ and $[S]_i$ represent the extracellular and intracellular concentrations
of ion $S$. Values for these constants are provided in Table \ref{table:constants}, while concentration
values for select ion species as well as their Nernst potentials are listed in
Table \ref{table:concentrations}. 

In summary, the combination of the diffusion and electrical forces is called the \emph{electrochemical gradient},
and can be thought of as a form of potential energy for a cell. In fact, the energy required for an action 
potential (discussed briefly in Chapter 1) is provided almost entirely from the electrochemical gradient
between the extracellular and intracellular regions of a neuron. 

\begin{table}
\begin{center}
  \caption{Various ionic concentrations, Nernst potentials, and resting
potentials for different cell types. Voltages are in $mV$ and concentrations in
$mM$. Adapted from \cite{keener-2009}.}
  \label{table:concentrations}
  \begin{tabular}{ | p{2.7cm} p{1.5cm} p{2.0cm} p{1.9cm} |}
    \hline
       & Squid Giant Axon & Frog \mbox{Sartorius} Muscle & Human Red
Blood Cell \\ \hline
      Intracellular concentrations & & & \\ \hline
      $Na^+$ & 50 & 13 & 19 \\
      $K^+$ & 397 & 138 & 136 \\
      $Cl^-$ & 40 & 3 & 78 \\ \hline
      
      Extracellular concentrations & & & \\ \hline
      $Na^+$ & 437 & 110 & 155 \\
      $K^+$ & 20 & 2.5 & 5 \\
      $Cl^-$ & 556 & 90 & 112 \\ \hline    
      
      Nernst \mbox{potentials} & & & \\ \hline
      $V_{Na}$ & +56 & +55 & +55 \\
      $V_{K}$ & -77 & -101 & -86 \\
      $V_{Cl}$ & -68 & -86 & -9 \\ \hline   
      
      Resting & & & \\
      potentials & -65 & -99 & -6 to -10 \\ \hline
  \end{tabular}
\end{center}
\end{table}

\subsection{The cell membrane}
Biological cells are divided into two main groups: prokaryotic cells which
include bacteria, and eukaryotic cells such as plant and animal cells.
Prokaryotic and eukaryotic cells are both enclosed in a selectively permeable
biological membrane known as the cell (or plasma) membrane. Eukaryotic
cells also have a
number of internal membranes enclosing their nucleus, mitochondria, and
several other organelles. Prokaryotes on the other hand do not have a nucleus or
any other membrane-bound organelles. 

The cell membrane of both prokaryotes and
eukaryotes, as well as the membranes of certain organelles all contain ion
channels. From this point on, however, the focus will be on the cell membrane
of eukaryotic cells. 

The cell membrane separates the intracellular contents from the
extracellular environment and is composed of a thin double layer of lipids about
7.5 $nM$ thick \cite{keener-2009}. A lipid is a simple molecule with one
hydrophilic head and one or more hydrophobic tails. Lipids group together in a
double layer to form the phospholipid bilayer of the cell membrane. Dispersed
throughout the cell's membrane are an an assortment of integral (also called
trans-membrane) proteins, as well as a variety peripheral proteins and other
molecules such as carbohydrates. The variety and importance of membrane proteins 
is evident by a recent study estimating that such proteins
account for $~30\%$ of genes \cite{capener-2002}. Figure
\ref{fig:cell_membrane} provides a view of a typical cell membrane. 

%%%% Cell membrane  %%%%
\begin{figure}[h!]
\begin{center}
 
\fbox{\includegraphics[width=5.8in]
{cell_membrane.eps}}
    \caption{Representation of the cell membrane showing the lipid bilayer,
integral and peripheral proteins, and carbohydrates. From Wikipedia \cite{wiki-cellmembrane}.}
    \label{fig:cell_membrane}
\end{center}
\end{figure}

An important characteristic of the cell membrane is that it is selectively
permeable, meaning that it allows certain molecules to pass though in either
direction in controlled quantities. For example, waste products are allowed to
leave the cell, and certain sugars, acids, and other molecules can enter.
Also affected by the membrane's selective permeability are ions, such as
$Na^+$ (sodium), $K^+$ (potassium), $Cl^-$ (chloride), and $Ca^{2+}$ (calcium).

The transportation of molecules across the membrane can either be passive
(requiring no net energy), or active (requires energy, such as ATP). An
important active transport mechanism is the
$Na^+/K^+-ATPase$, which is also called the sodium-potassium pump for the
sake of readability. This pump is a trans-membrane protein that uses the
energy stored in ATP to pump $Na^+$ out of the cell and $K^+$ into the cell
\cite{keener-2009}. In addition to pumps, integral proteins called exchangers
also facilitate active transport across the cell membrane. The $Na^+-Ca^{2+}$
exchanger uses the energy gained from allowing $Na^+$ entry into the cell to
remove $Ca^{2+}$. One may question why allowing $Na^+$ into the cell provides
energy, while removing $Ca^{2+}$ requires energy. The reason this is so is due
to the electrochemical gradients arising from the ionic concentration levels
that the cell maintains through active transport. 

At the cell's resting potential, the concentration of $Na^+$ is much higher
outside the cell than inside the cell, so $Na^+$ wants to enter because the
electrochemical gradient is driving it into the cell.
$Ca^{2+}$ also a higher extracellular
concentration, so it takes energy to remove it from the cell because it must be
moved up, or \emph{against} its electrochemical gradient.


The last transport mechanism that deserves mention are ion channels,
trans-membrane proteins which facilitate the diffusion of select ion species
down their electrochemical gradient. Ion channels are therefore
referred to as passive transport mechanisms. As a simplification, an individual channel 
can be thought of as either open, or closed. 
In the open configuration, channels allow ions to pass
through at speeds of $~10^8$ ions per second, near diffusion speed
\cite{capener-2002}.

From the perspective of the whole cell, the large number ion channels together
can be said to affect the \emph{permeability} of the membrane to a certain ion
species.  

In summary, the cell membrane uses ion channels, pumps, and exchangers to
maintain a certain concentration gradient for a number of ion species. The
membrane's permeability to an ion species depends on the types of ion
channels in that cell type, as well as their current states (open, or closed).


\subsection{Membrane potential}
Given various intracellular and extracellular ionic concentrations, as well as
the membrane's relative permeability to these ions, one can calculate the
resting membrane potential of a cell. The resting membrane potential  First,
some basic concepts in
electricity need to be refreshed.

Quantity of charge is measured in coulombs, and is abbreviated $C$. Current, ($I$),
represents the movement of charge and has units coulombs per second, or amperes
($A$). Conductance, ($G$), is the inverse of resistance ($R$), and measures the
ease of current flow between two points. The unit of conductance is siemens
($S$), while the unit of resistance is ohms ($\Omega$). Voltage, $V$, measures
the potential difference between two points. Note that the terms ``potential'' and
``voltage'' can be used interchangeably.  

Ohm's law defines the relationship between voltage, current, and resistance:
\begin{equation}
 V = I R.
\end{equation}

It may also be written in terms of conductance:
\begin{equation}
 V = \slfrac{I}{G}.
\end{equation}


Now let us apply Ohm's law to a biological cell. 
Considering the simplification that the membrane is permeable only 
to ion species $S$, the membrane's
conductance, $G$, is directly related to how easily $S$ passes though its ion
channels. This is called channel conductance, and is a function of the
opening and closing (referred to as gating from now on) of the channels as well
as their type. $I_{S}$ is the current through the membrane due to the
movement of species $S$. In this case, the voltage, $V$, is referred to as
the \emph{driving force}. The driving force is different for each different ion
species in the cell - $Na^+$ ions, for example, experience a force greater in
magnitude and opposite in direction compared with $K^+$ ions. 

It is worth
noting that it is difficult to measure the conductance of a cell. Permeability, 
on the other hand, is easy to measure. Permeability and conductance are related but 
not equivalent; however, they are treated as such in electrophysiology.
 	 

To summarize: the current from ion $S$ is equal to the driving force of $S$
times the membrane conductance to $S$. Formally:
\begin{equation}
 I_{S} = G_{S} (DrivingForce_S).
\end{equation}

Driving force is determined by the membrane voltage, $V_m$, and $E_S$, the
equilibrium potential (calculated using Equation \ref{eq:nernst}) for ion $S$:
\begin{equation}
 DrivingForce_S = (V_m - E_S).
\end{equation}

Substituting:
\begin{equation}
 I_{S} = G_{S} (V_m - E_S).
 \label{eq:current}
\end{equation}
Equation \ref{eq:current} is used by the software described in Chapter 3 to calculate the
ionic currents given a known membrane voltage and simulated channel
conductances.

To calculate the resting potential of a cell permeable to more than one ion type,
we need to take into account the
concentrations of the various ion species, as well as the conductance of their
respective channels. The equation that does this is the Goldman-Hodgkin-Katz
(GHK) voltage equation:

\begin{equation}
 E_{m}=\frac{RT}{F} \; ln \frac{P_K[K]_e + P_{Na}[Na]_e + P_{Cl}[Cl]_i}
  {P_K[K]_i + P_{Na}[Na]_i + P_{Cl}[Cl]_e},
\end{equation}

where $E_{m}$ is the resting membrane potential and $P_{S}$ represents the
relative permeability of the membrane to ion species $S$. 

It is worth stating that unlike the Nernst equation which is universally true,
the GHK equation depends on the assumption of a constant electric field.


\begin{table}
\begin{center}
  \caption{Physical constants. From \cite{hille-2001}.}
  \label{table:constants}
  \begin{tabular}{ | l l l l |}
    \hline
      Faraday's constant & $F$ & $9.6485 \times 10^4$ & $Cmol^{-1}$ \\  
      Absolute temperature & $T(K)$ & $273.15$ & $^\circ Celcius$ \\
      Gas constant & $R$ & $8.3145$ & $VCmol^{-1}K^{-1}$ \\
    \hline
  \end{tabular}
\end{center}
\end{table}


%\subsection{Membrane ion channels}

%\subsection{Electrophysiology}
%\label{section:electrophysiology}

%%%% Patch clamp  %%%%
%\begin{figure}[h!]
%\begin{center}
%\includegraphics[width=3in]
%{patch_modes.eps}
%    \caption{Various patch clamp techniques. From Wikipedia \cite{}}
%    \label{fig:patch_clamp}
%\end{center}
%\end{figure}

\section{Modeling ion channels}
The modeling of ion channels dates back over 100 years
\cite{levitt-1999}. As
is with many mathematical models, the motivations driving the development of ion
channel models include:

\begin{itemize}
 \item striving to infer fundamental knowledge of the underlying physical
processes,
 \item matching experimental data to the model's predictions,
 \item using the model to predict behaviors which are difficult to observe, and
 \item developing common terminology and knowledge among researchers.
\end{itemize}

The interpretation of ion channels models during the pre-molecular era provided
the primary source of information about channel structure \cite{levitt-1999}.
However, current technology such as molecular biology and x-ray diffusion are
beginning to prove atomic-resolution detail of bacterial as well as some
mammalian channels, though the fraction of channels with known structure
remains low \cite{levitt-1999, capener-2002}. 

Hodgkin and Huxley's landmark model describing the ionic mechanisms behind
the action potential of a squid giant axon provided the foundation for many
modern practices in electrophysiology and ion channel modeling
\cite{hodgkin-1952, keener-2009}. Their model describes the voltage-dependent
and time-dependent behavior of $Na^+$ and $K^+$ conductances using a coupled set
of ordinary differential equations. While the Hodgkin-Huxley
model is still used today due to its simplicity and low number of parameters, it
does exhibit several shortcomings. These include, for example, a lack
of connectivity between activating and inactivating gates in the $Na^+$ channel,
as well as the premise that the inactivation gate can only close after the
activation
gate opens \cite{gurkiewicz-2007}.

\subsection{Markov models}
Markov models for ion channels are an extension of the Hodgkin-Huxley formalism, 
and are good for modeling single channel data, gating currents, and drug
interactions. Unlike Hodgkin-Huxley models, Markov models show the state
dependence of activation and inactivation \cite{fink-2009}. 

\subsubsection{Theory}
Ion channels can be modeled as continuous time Markov processes (see Section
\ref{section:markov}), where a simplification of the channels functional
physical shape is represented by states in the Markov chain. A model can have
several open, closed, deactivated, and inactivated states, as well as states
which represent the binding of ligands. The rates connecting the various
states are time independent kinetic rates which can be constant, or dependent on
the membrane voltage, or intra- and extra-cellular ionic concentrations.

\subsubsection{Criticism}
Markov models are not without critics. In a one page perspective
titled \emph{Are rate constants constant?} Jones challenges the community to
examine the benefits of time-dependent rate constants over Markov models
\cite{jones-2006}. Jones cites the findings by Uebachs \cite{uebachs-2006} who
reports that the rate of recovery from inactivation in a type of
$Ca^{+2}$ channel depends on the length of depolarization, which goes against
the basic principle of chemical kinetics that rate constants remain constant in
a constant environment \cite{jones-2006}.

Although Markov models are memoryless by definition, the time-dependent behavior
described
by Uebachs can be reproduced by having a long chain of open or closed states
\cite{jones-2006, liebovitch-1989}. Nevertheless, Jones concludes that
time-dependent rate constants as used by fractal channel models, would provide a
more explicit memory than multistate Markov models.  

\subsection{Fractal models}
Liebovitch has been the main proponent of fractal ion channel models for over
15 years \cite{liebovitch-1989, lowen-1999, liebovitch-2001}. One of his main arguments questioning the use
of Markov models is founded in the body of evidence depicting ion channel
proteins as complex structures with a continuum of states, exhibiting large and
small motions over varying time scales. Fractal models as proposed by Liebovitch
contain a large number of conformational states with kinetic rates connecting
them. The rates, however, are linked and not independent as in Markov models.
This has the effect of giving the channel ``memory'' \cite{liebovitch-1989}. 

An additional concern raised by Liebovitch is the use of exponential rate
constants. He provides a friendly reminder to his colleagues
that the fitting of experimental data to sums of exponentials is famously
ill-conditioned \cite{liebovitch-1989, acton-1970}. As confirmation he provides
a quote from a numerical methods text which states that those who try to
determine the parameters of such equations from experimental data ``must be
spanked or counseled. At the very least, keep them from obstructing Progress and
the computer!'' \cite{acton-1970}. 

In their defense, other researchers have argued that Markov models with
exponential rate constants and a small number of states fit the experimental
data better than fractal models \cite{mcmanus-1988}. Nevertheless, fractal
models continue as the number one challenger to Markov models.

\section{Math Background}

\subsection{Continuous time Markov process}
\label{section:markov}
The following information is inspired from David Anderson's notes on
biological stochastic processes \cite{anderson-2011}. 

Continuous time Markov processes are commonly used to simulate stochasticity in
a variety of fields including biology, chemistry, physics, and economics. As
the name implies, a continuous time Markov process is similar to a discrete
time Markov chain, with the exception that transitions between states can occur
at any time, with a exponentially distributed probability. Both discrete and
continuous time Markov processes satisfy the Markov property which says that
the probability distribution of future states depends only on the current state,
and not on the preceding events. This ``memoryless'' behavior is what gives
Markov processes their stochastic behavior. The property is named after the
Russian mathematician Andrey Markov. 

Continuous time Markov chains are easily described using weighted directed
graphs, where nodes represent the states, and the edges are labeled with the
\emph{rates} of transitioning between states. Figure
\ref{fig:three_state_markov} shows a three state Markov chain in which state $a$
transitions to state $b$ with rate $k_1$, and $c$ transitions to $b$ with a rate
$k_4$. Note that unlike the discrete case, the transition rates are not
equivalent to probabilities and thus are not restricted to values $[0-1]$.

%%%% Markov chain  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=3in]
{three_state_markov.eps}
    \caption{A three state Markov chain with states a, b, and c, and
transition rates $k_1 - k_4$.}
    \label{fig:three_state_markov}
\end{center}
\end{figure}


To represent the connections between $N$ states in a continuous time Markov
chain, we can use a weighted adjacency
matrix, $D$, where $D$ is size $N \times N$, and the non-diagonal entries,
$D_{ij}$ represent the edge weight, or \emph{transition rates} from state $j$
to state $i$. Self-connections are not allowed, i.e., state $b$ cannot define a
rate to itself. Therefore, the diagonal entries of the matrix must be zero.

The adjacency matrix representing the continuous time Markov chain given in
Figure \ref{fig:three_state_markov} is:

\begin{equation}
D = 
\begin{bmatrix}
0  & k_3 & 0 \\
k_1 & 0 & k_4 \\
0 & k_2 & 0 \end{bmatrix}
\end{equation}

The adjacency matrix describes the connections of a Markov chain. To describe
the rates in which we move between states, we use the transition matrix
(also called the infinitesimal
generator matrix), $A$, which is the adjacency matrix with special entries in
the diagonal such that each column sums to zero \cite{anderson-2011}. 
Thus the diagonal entries must
be equal to the negative sum of all the rates for that column. For the
three state example, we get the following transition matrix:

\begin{equation}
A = 
\begin{bmatrix}
-k_1  & k_3 & 0 \\
k_1 & -(k_2+k_3) & k_4 \\
0 & k_2 & -k_4 
\end{bmatrix} 
\end{equation}

Now we can describe the evolution of the state probabilities by deriving a set
of ordinary differential equations, called the Kolmogorov forward equations \cite{anderson-2011}. 
In the case of natural sciences, they are known as the master equation, or
chemical master equation. The master equation has the form:

\begin{equation}
\label{eq:master}
\frac{d\vec{P}}{dt} = A\vec{P},
\end{equation}

where $\frac{d\vec{P}}{dt}$ represents the change of the
state probability vector $\vec{P}$, and $A$ is the transition matrix. The
master equation in its expanded form for the three state example is:

\begin{equation}
\begin{bmatrix}
 \slfrac{dP_a}{dt} \\
 \slfrac{dP_b}{dt} \\
 \slfrac{dP_c}{dt} 
\end{bmatrix} = 
\begin{bmatrix}
-k_1  & k_3 & 0 \\
k_1 & -(k_2+k_3) & k_4 \\
0 & k_2 & -k_4 
\end{bmatrix} 
\begin{bmatrix}
 P_a \\
 P_b \\
 P_c 
\end{bmatrix} = 
\begin{bmatrix}
 -k_1 P_a + k_3 P_b \\
 k_1 P_a - (k_2+k_3)P_b + k_4P_c \\
 k_2P_b - k_4P_c 
\end{bmatrix}
\end{equation}

We now have three differential equations to describe the evolution of our
three state system. However, this system is over specified. Instead of
explicitly calculating the evolution of the last state (state $c$ in our
example), we infer it via the conservation of probabilities like so:

\begin{equation}
 P_N = 1- \displaystyle\sum_{i=1}^{N-1} P_i.
\end{equation}

In other words, the probability of model being in the last state, plus the sum
of the other state probabilities must be equal to
one.

Using the conservation of probabilities, our system can be re-written as:

\begin{equation}
\begin{bmatrix}
 \slfrac{dP_a}{dt} \\
 \slfrac{dP_b}{dt} \\
 P 
\end{bmatrix} = 
\begin{bmatrix}
-k_1  & k_3 & 0 \\
k_1 & -(k_2+k_3) & k_4 \\
1 & 1 & 1 
\end{bmatrix} 
\begin{bmatrix}
 P_a \\
 P_b \\
 P_c 
\end{bmatrix} = 
\begin{bmatrix}
 -k_1 P_a + k_3 P_b \\
 k_1 P_a - (k_2+k_3)P_b + k_4P_c \\
 P_a + P_b + P_c 
\end{bmatrix}.
\end{equation}

If we wish to solve our system of ordinary differential equations as an initial
value problem, we need initial conditions. When simulating the current through
ion channels, it is often desired to start
solving the system while it is at steady state, that is, when all of the state
probability derivatives are equal to zero:

\begin{equation}
\begin{bmatrix}
 \slfrac{dP_a}{dt} \\
 \slfrac{dP_b}{dt} \\
 P 
\end{bmatrix} = 
\begin{bmatrix}
-k_1  & k_3 & 0 \\
k_1 & -(k_2+k_3) & k_4 \\
1 & 1 & 1 
\end{bmatrix} 
\begin{bmatrix}
 P_a \\
 P_b \\
 P_c 
\end{bmatrix} = 
\begin{bmatrix}
 0 \\
 0 \\
 1
\end{bmatrix}, or
\end{equation}

\begin{equation}
 \vec{P^\prime} = A \vec{P_{ss}} = \vec{P}.
\end{equation}

Thus,

\begin{equation}
\label{eq:ss}
 \vec{P_{ss}} = \slfrac{A}{\vec{P}}.
\end{equation}

Equation \ref{eq:ss} is used in the implementation of our software
(see Chapter 3) to set the state probabilities to steady state before
beginning integration. During simulation we supply the derivatives of the state
probabilities (calculated using Equation \ref{eq:master}) to the solver during each
integration step.

%\clearpage
\section{Existing ion channel simulators}
We provide a comparison of several Markov model based ion channel simulators.
Table \ref{table:simulators} compares the basic features of each simulator,
including our solution, ModFossa.

\begin{table}
\scriptsize
\begin{center}
  \caption{Feature comparison of several ion channel simulators which use
Markov models.}
  \label{table:simulators}
  \begin{tabular}{ | p{1.9cm} | p{0.8cm} |p{0.4cm} |p{1.2cm} |p{1.2cm} |p{0.8cm}
|p{1.8cm} |p{0.6cm} |p{1.5cm} |} \hline
    Name & Model & GUI & Scriptable & Plotting & Plot \mbox{format} & Simulation
type & OS & License \\ \hline
    IonChannelLab & Markov & yes & no & yes, \mbox{interactive} & bitmap &
    numerical integration, Monte Carlo, Q-Matrix & Win. & Free, no source
    available \\ \hline
    ChannelLab & Markov & yes & no & yes & bitmap & numerical
    integration, Monte Carlo & Win. & Proprietary \\ \hline
    QUB & HMM & yes & yes & yes, \mbox{interactive} & bitmap & numerical \mbox{integration} & Win.,
    Mac, Linux & GPL \\ \hline
    ModFossa & Markov & no & yes, from Python & yes & vector & numerical
    \mbox{integration} & Linux & GPL \\ \hline
  \end{tabular}
\end{center}
\end{table}



\subsection{IonChannelLab}
IonChannelLab is a standalone Windows application targeted towards the
non-programmer, and is therefore focused on ease of use \cite{castillo-2010}.
It is available free of charge; however, source code is not available.

\subsubsection{General observations}
The software provides a nice looking interface which is uncluttered, intuitive,
and easy to learn; however, the lack of online help 
(especially context-sensitive help) proves somewhat disappointing. 
Several example projects are available, including Hodgkin-Huxely type models, as well
implementations of the author's own research on chloride channels.

IonChannelLab allows users to save projects in an XML format for later use.
Lastly, the software is not scriptable. If it were, IonChannelLab would no doubt
be more attractive to advanced users.

\subsubsection{Model definition}
IonChannel lab facilitates easy Markov model creation by allowing the user to 
drag and drop states into arbitrary locations on a 2D grid. A transition is created
by connecting two states using the mouse. Each transition specifies two rate constants, 
one for the rate from state A to state B, and another from B to A. Rate constants are 
specified by a unique string name, and are defined using the function editor 
window (Figures \ref{fig:icl_function_exp} and \ref{fig:icl_function_custom}).
The function editor allows the user to specify the parameters of several 
available rate constant types, including exponential
voltage; exponential voltage with temperature dependence; double exponential
with temperature dependence; ligand gated; exponential ligand gated with
temperature dependence; Hodgin-Huxely types A, B, and C; constant; custom
voltage; and custom ligand. The equation used to calculate each rate constant
type is clearly shown in the function editor window, a useful feature not found  
in other simulators reviewed here.

IonChannelLab's model definition screen also provides the user with two choices
for the permeation model, Ohm's law and Goldman-Hodgin-Katz.

%%%% Icl Model  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=6in]
{icl_model.png}
    \caption{IonChannelLab model definition screen. The main area contains the
Markov model with 12 states connected by various rate constants.
The rate constants and their equations are shown in the sidebar on the left.}
    \label{fig:icl_model}
\end{center}
\end{figure}

%%%% Icl function exp  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=4in]
{icl_function_exp.png}
    \caption{IonChannelLab function editor window showing an exponentially
voltage dependent rate constant to which the user has assigned the name
``a1''. The equation used to calculate the rate is shown in easy-to-read text.}
    \label{fig:icl_function_exp}
\end{center}
\end{figure}

%%%% Icl function custom  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=4in]
{icl_function_custom.png}
    \caption{IonChannelLab function editor window showing a custom rate
constant. The user has defined this particular rate constant as exponentially
voltage dependent.}
    \label{fig:icl_function_custom}
\end{center}
\end{figure}

\subsubsection{Stimulus definition}
Voltage and concentration protocols are defined in the stimulus window as
shown in Figure \ref{fig:icl_stimulus}. Constant, stepped, and ramped stages
are supported for both voltage and concentration protocols. 
A visualization of the protocol is provided in the stimulus
window, and is updated immediately after a change in any of the fields.

%%%% Icl stimulus  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=6in]
{icl_stimulus.png}
    \caption{IonChannelLab stimulus editor window showing a visualization of the
voltage protocol. As shown in the top of the window, this experiment also defines 
internal and external chloride concentration protocols.}
    \label{fig:icl_stimulus}
\end{center}
\end{figure}


\subsubsection{Experiment definition}
The experiment definition window provides settings and visualizations specific 
to the experiment. In creating a new experiment using IonChannelLab,
the user chooses from the following experiment types: current vs time, 
conductance vs time, gating current vs time, math expression vs time,
steady state current vs X, steady state conductance vs X, steady state 
expressions vs X, and time constants vs X.

The software gives the user the option of selecting one of several
different solvers, including q-matrix, integration via Runge-Kutta, integration 
via Gear's BDF, and single channel Monte Carlo. Implementation of the
solvers is provided by DotNumerics, a numerical library written using Microsoft's
.NET framework. The initial conditions of the simulation can be 
calculated from the holding values, or set manually by the user.
Additional experiment settings include the number of channels to  simulate, 
temperature, and an optional noise.

\subsubsection{Plotting and analysis}
All plots are interactive and allow for the user to pan and zoom.
If the user wishes to save a plot, he can either save the data in CSV format, 
or copy the image to the clipboard. There is no option
to export plots in a vector graphics format. 

\subsection{ChannelLab}
ChannelLab is a powerful, but dated commercial single ion channel modeling program developed
over the last 11 years by Synaptosoft, Inc. \cite{channellab}. It is available only for 
Microsoft Windows.

\subsubsection{General observations} 
Unlike the other two simulators reviewed here, ChannelLab
is not free; an individual license costs \$250. While the 
software is quite old (the last version was released in 2003), 
it contains many useful features, such as curve fitting 
and analysis tools, and data conversion utilities. 
Over 100 published channel models are included as examples.

Nevertheless, the interface does appear dated, especially when adding transitions between
states. Additionally, ChannelLab lacks any scripting support.

\subsubsection{Model definition}
ChannelLab supports up to 20 states, with each state having an adjustable 
gating parameter. To add connections between states, the user can either use
the 2D grid layout (Figure \ref{fig:cl_grid}), or the 3D table layout 
(Figure \ref{fig:cl_table}). Both layouts are decidedly less intuitive and more 
constraining than the other simulators mentioned here. 

Rate constants support dependence on up to six drugs or voltages, although only
the user is limited to exponential equations.

%%%% Cl grid  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5in]
{cl_grid.png}
    \caption{ChannelLab 2D grid layout of a four state model connected with
    six rate constants.}
    \label{fig:cl_grid}
\end{center}
\end{figure}

%%%% Cl table  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5in]
{cl_table.png}
    \caption{ChannelLab 3D table layout of the same model shown in Figure \ref{fig:cl_grid}.}
    \label{fig:cl_table}
\end{center}
\end{figure}

\subsubsection{Stimulus definition}
The software provides a powerful stimulus definition window, as shown in Figure \ref{fig:cl_stimulus}.
Concentration and voltage protocols are defined using automatic step functions, manual step or ramp
functions, exponential functions, or stimulus trains. The protocols are plotted in the stimulus
definition window using either linear or log scaling in order to provide immediate feedback to the user.
   
Channel lab also allows for the import stimulus and voltage waveforms from ascii text files.

%%%% Cl stimulus  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=6in]
{cl_stimulus.png}
    \caption{The ChannelLab stimulus window allows the user to create complex voltage and concentration
    protocols while providing immediate visual feedback.}
    \label{fig:cl_stimulus}
\end{center}
\end{figure}

\subsubsection{Experiment definition}
Three solvers are supported by ChannelLab: Runge-Kutta 4, Runge-Kutta 5, and Monte Carlo.
As expected, the software allows the user to select the initial conditions as either the steady
state probabilities, or user-provided values.


\subsubsection{Plotting and analysis}
The most attractive features of ChannelLab are the interactive analysis tools, which
include maximum likelihood and least squares waveform fitting. A screenshot
of the analysis window is provided in Figure \ref{fig:cl_analysis}. The plots in the 
analysis window are interactive and can be exported to a variety of file formats, but cannot
be saved as vector graphics. 


%%%% Cl analysis  %%%%
\begin{figure}[h!]
\begin{center}
\includegraphics[width=6in]
{cl_analysis.png}
    \caption{ChannelLab analysis window. Note the large number of adjustable analysis settings to the right.}
    \label{fig:cl_analysis}
\end{center}
\end{figure}


% \subsection{QUB}
% \subsubsection{General observations}
% 
% Window layouts a little bit confusing at first. Realtime feedback. 
% 
% States are color coded depending of their conductance. 
% \cite{qub}
% 
% \subsubsection{Model definition}
% \subsubsection{Stimulus definition}
% \subsubsection{Experiment definition}
% \subsubsection{Plotting and analysis}
% 
% 
% %%%% Qub Main %%%%
% \begin{figure}[h!]
% \begin{center}
% \includegraphics[width=6in]
% {qub_main.png}
%     \caption{Qub express main window. The bottom panel shows the model area which consists of states (boxes) 
%     and connections (arrows). Above the model area is the table planel where all aspects of the experiment are 
%     defined, including states, rates, and stimuli. The middle panel shows the simulated data from the current 
%     model. The top panel serves various purposes, including curve fitting, and administration.  
%     }
%     \label{fig:qub_main}
% \end{center}
% \end{figure}
% 
% Summary
% \begin{itemize}
%  \item Windows, Mac, Linux
%  \item Free
%  \item Source available on Git
%  \item Uses Python, boost++
%  \item Currently maintained, updates this month
%  \item Hidden Markov model
%  \item Statistical clustering (K-Means)
%  \item Separate application for data acquisition
%  \item Can record actions as Python script
%  \item Host plugins
% \end{itemize}



