%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Trim Size: 9.75in x 6.5in
%% Text Area: 8in (include Runningheads) x 5in
%% ws-jbcb.tex   :   24-5-04
%% Tex file to use with ws-jbcb.cls written in Latex2E. 
%% The content, structure, format and layout of this style file is the 
%% property of World Scientific Publishing Co. Pte. Ltd. 
%% Copyright 1995, 2002 by World Scientific Publishing Co. 
%% All rights are reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

\documentclass{ws-jbcb}
\usepackage{url} %added fch
%\usepackage{amsmath} % This command was added to ws-jbcb.cls by Ferneyhough
\usepackage{graphicx} % added by fch
\usepackage{epstopdf} % added fch
\usepackage{verbatim} % added fch

% Right slash in math mode for fractions.
\newcommand{\slfrac}[2]{\left.#1\middle/#2\right.}

% Custom colors
\usepackage{color}
\definecolor{deepblue}{rgb}{0,0,0.5}
\definecolor{deepred}{rgb}{0.6,0,0}
\definecolor{deepgreen}{rgb}{0,0.5,0}
\definecolor{light-gray}{gray}{0.95}

\usepackage{setspace}
\usepackage{listings}
\usepackage{url}

% Python style for external listing
\newcommand\pythonstyleA{
\lstset{
language=Python,
basicstyle=\scriptsize\ttfamily,
otherkeywords={self},
keywordstyle=\ttb\color{deepblue}\scriptsize\ttfamily,
emph={MyClass,__init__},
commentstyle=\ttb\color{deepred}\scriptsize\ttfamily,
stringstyle=\color{deepgreen}\scriptsize\ttfamily,
frame=single,
showstringspaces=false
}}

% Python style for in source listing
\newcommand\pythonstyleB{
\lstset{
language=Python,
backgroundcolor=\color{light-gray},
basicstyle=\ttm,
otherkeywords={self, True, False},
keywordstyle=\ttb\color{deepblue},
emph={MyClass,__init__},
commentstyle=\ttb\color{deepred},
stringstyle=\color{deepgreen},
showstringspaces=false
}}

% Python environment
\lstnewenvironment{python}[1][]
{
\pythonstyleB
\begin{singlespace}
\lstset{#1}
\end{singlespace}
}
{}

% Python for external files
\newcommand\pythonexternal[1]{{
\pythonstyleA
\begin{singlespace}
\lstinputlisting[float=h]{#1}
\end{singlespace}}
}

% Python for inline
\newcommand\pythoninline[1]{{\pythonstyleB\lstinline!#1!}}


\begin{document}

\markboth{Ferneyhough, Thibealut, Dascalu, Harris}
%{Modfossa: A Python Library for Ion Channel Modeling}
{Modfossa: A Library for Modeling Ion channels using Python}

%%%%%%%%%%%%%%%%%%%%% Publisher's Area please ignore %%%%%%%%%%%%%%%
%
\catchline{}{}{}{}{}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{ModFossa: A Library for Modeling Ion channels using Python}
% \footnote{For the title, try not to 
% use more than 3 lines. Typeset the title in 10~pt 
% Times roman, uppercase and boldface.} }

\author{Gareth B. Ferneyhough, Corey M. Thibealut, Sergiu M. Dascalu, and Frederick C. Harris, Jr.}
\address{Department of Computer Science and Engineering,\\
University of Nevada, Reno,\\
Reno, Nevada 89557,
USA\\
Fred.Harris@cse.unr.edu}

\maketitle

\begin{history}
\received{(Day Month Year)}
\revised{(Day Month Year)}
\accepted{(Day Month Year)}
%\comby{(xxxxxxxxxx)}
\end{history}

\begin{abstract}
The creation and simulation of ion channel models using continuous-time Markov
processes is a powerful and well-used tool in the field of electrophysiology
and ion channel research. While several software packages exist for the purpose
of ion channel modeling, none are available as a Python library. In an attempt
to provide an easy-to-use, yet powerful Markov model-based ion channel
simulator, we have developed ModFossa, a Python library supporting easy model
creation and stimulus definition, complete with a fast numerical solver, and
attractive vector graphics plotting. 
\end{abstract}

\keywords{Ion channel modeling; Markov models; Python.}

\section{Introduction}
Ion channels are trans-membrane proteins that control the passive flow of ions
into, and out of a biological cell. They exist in all cell types and play
crucial roles in cellular sensing and communication, maintenance of the
membrane potential, and regulation of cell volume.

For example, voltage-gated sodium, potassium, and calcium channels, among
others, work together to shape the automatically recurring cardiac action
potential that give
our hearts their sure, and usually steady, beat. Philosophy aside, the very
thoughts passing through the consciousness of a scientist as he reads this page
exist as nothing more than electro-chemical signals in the brain. Specifically,
voltage-gated sodium and potassium channels integrated in the cell membrane
of all neurons allow electrical signals to propagate rapidly down the axon in
the form of an action potential. When the action potential reaches the end of
the axon, voltage-gated calcium channels regulate the release of special
chemical messengers called neurotransmitters into the synapse,the structure
that connects neurons to one another.
As the neurotransmitters diffuse across the synapse they bind to ligand-gated
ion channels on the post-synaptic neuron. Depending on the type of ion channel
and the type of neurotransmitter, the channels may open or close, thereby
altering the membrane voltage. If
the post-synaptic neuron receives enough excitatory stimulus from the
many presynaptic neurons in which it is connected, its membrane potential will
depolarize to a point and it too will experience an action potential,
thus passing the signal along to other neurons.  

As further evidence of their vital importance, ion channels are specifically
targeted by several naturally occurring toxins and poisons, including venom
from spiders, snakes, scorpions, and fish. It comes as no
surprise that the development of new types of drugs which alter the behavior
of certain channels is a major motivator behind ion channel research. As in most
scientific fields, researchers often develop mathematical models to assist in
their understanding of the results. Ion channels are no exception; several
classes of ion channel models exist that allow for the study of ion channel
behavior under varying conditions. One such generalization involves modeling
channel
gating using continuous time Markov processes, and solving the underlying
ordinary differential equations that arise. A number of commercial as
well as open-source software programs exist
for this purpose, but it is also common for scientists to implement their own
solution using the differential equations solvers in Matlab, for example. 

In order to provide an elegant and simple solution for the creation and
simulation of Markov model-based ion channels, we present ModFossa, a simulator
written in C++ with an easy-to-use Python interface. Several published models
of ion channel gating have been reproduced using our simulator; these results,
along with performance analysis are presented in the following chapters.
In Section 2 we present several mathematical, biological, and electrical
concepts related to ion channel modeling. In Section 3, we present the design
and implementation of
ModFossa. Results are provided in Section 4. Finally, we conclude with a
discussion of applications and future work in Section 5. 

\section{Background}

\subsection{Biology background}

\paragraph{Salts, ions, and the cellular solutions:}
Ions, charged atoms, arise from a variety of natural processes, including the dissolution of salts
in water. For example, 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. 

To see the importance of ions in the role of cellular function, imagine a bath containing two 
solutions separated by a membrane which is perfectly permeable to $K^+$. 
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. 
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$. 

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. 

\paragraph{Cell membrane:}
The cell membrane separates the intracellular and extracellular regions of a cell
and is composed of a thin double layer of lipids about
7.5 $nM$ thick.\cite{keener-2009}  An important characteristic of the cell membrane is that it is selectively
permeable, which is facilitated through several mechanisms including ion channels -
trans-membrane proteins which facilitate the passive diffusion of select ion species
down their electrochemical gradient. 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.  

\paragraph{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. First, some basic concepts in 
electricity need to be refreshed.

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

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. 

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 Section 3 to calculate the
ionic currents given a known membrane voltage and simulated channel
conductances.

\subsection{Modeling ion channels}

\begin{comment}
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}

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

\subsection{Math Background}

\subsubsection{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. 

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}
\begin{center}
\includegraphics[width=1.6in]
{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
$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 Section 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.

\end{comment}
%---------------------------------------------

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:
(1) striving to infer fundamental knowledge of the underlying physical
processes,
(2) matching experimental data to the model's predictions,
(3) using the model to predict behaviors which are difficult to observe, and
(4) developing common terminology and knowledge among researchers.


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{capener-2002,levitt-1999}

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}

{\bf Markov models:}
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. 
Continuous time Markov chains are easily described using weighted directed
graphs, where vertices represent the states, and the edges are labeled with the
\emph{rates} of transitioning between states. 

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}

{\bf Theory: }
Ion channels can be modeled as continuous time Markov processes, 
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.

{\bf 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.  

{\bf Fractal Models:}
Liebovitch has been the main proponent of fractal ion channel models for over
15 years.\cite{liebovitch-1989,liebovitch-2001,lowen-1999}  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 {\em 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{acton-1970,liebovitch-1989}  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.


%----------------------------------------------



\begin{comment}
\paragraph{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.

%\smallskip \noindent \textbf{General observations}\\
The software provides a nice looking interface which is uncluttered, intuitive,
and easy to learn; however, it lacks online and context-sensitive help.
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.

\begin{comment}
\smallskip \noindent \textbf{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.
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.

\smallskip \noindent \textbf{Stimulus definition}\\
Voltage and concentration protocols are defined in the stimulus window.
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.


\smallskip \noindent \textbf{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.

\smallskip \noindent \textbf{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. 



\paragraph{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.

%\smallskip \noindent \textbf{General observations} \\
Unlike the other two simulators reviewed here, ChannelLab
is not free. 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.


\smallskip \noindent \textbf{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
a 2D grid layout, or a 3D table layout. 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.

\smallskip \noindent \textbf{Stimulus definition}\\
The software provides a powerful stimulus definition window.
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.

\smallskip \noindent \textbf{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.


\smallskip \noindent \textbf{Plotting and analysis}\\
The most attractive features of ChannelLab are the interactive analysis tools, which
include maximum likelihood and least squares waveform fitting. 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. 
\end{comment}


\section{ModFossa}

\subsection{Overview}

In this section we present ModFossa, a fast and easy-to-use Python library for creating and
simulating ion channel kinetics using continuous time Markov processes. A
brief discussion of ModFossa's features, components, are provided here.  For
more details regarding the design of ModFossa the reader is referred to Ref.~\refcite{garreth-thesis} and for the source
code the reader is referred to Ref.~\refcite{garreth-source}.

For those interested, the name ModFossa was formed from two Latin words ---
mod, and fossa, which can be interpreted as \emph{open channel} when written
together.  
Before we begin the design overview,
Table \ref{table:simulators} compares the basic features of several simulators,
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.0cm} |p{1.2cm} 
|p{1.8cm} |p{0.6cm} |p{1.4cm} |} \hline
    Name & Model & GUI & Script & Plotting, \mbox{format} & Simulation
type & OS & License \\ \hline
    IonChannelLab & Markov & yes & no &  \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 & \mbox{interactive} bitmap & numerical \mbox{integration} & Win.,
    Mac, Linux & GPL \\ \hline
    ModFossa & Markov & no & yes,  Python &  yes, vector & numerical
    \mbox{integration} & Linux & GPL \\ \hline
  \end{tabular}
\end{center}
\end{table}


\subsection{Features}
%\paragraph{Python interface:}
{\bf Python interface:}
ModFossa provides easy, yet versatile and powerful ion channel modeling through
its novel Python interface. Python is well-used in the scientific and academic
community; however, we are unaware of any Python libraries which facilitate
the creation and simulation of ion channel kinetics. The popularity of Python
as a research tool stems from its simple syntax and extensive standard library,
and the availability of many high-quality open source numeric and scientific
libraries. Additional qualities of Python include easy third party library
installation; strong, dynamic typing; and an interactive interpreter.
ModFossa can also be used directly as a C++ library, allowing for more
flexibility.

%\paragraph{Easy model creation:}
%\noindent
{\bf Easy model creation:}
ModFossa allows the user to define an ion channel's Markov model using states,
rates, and connections. States and rates are referred to by their user-supplied
names. This allows the user to create the model in any order. For example,
connections can be defined before the definition of the states and rate
constants.

%\smallskip \noindent \textbf{States} \\
States represent the various channel states as defined
by a continuous time Markov chain. All states are referred to by unique name
which is assigned by the user during
creation. States are created as either non-conducting or conducting.
Additionally, the gating of conducting states can by specified during
construction. The gating parameter is used to simulate states that are
partially conducting. For example, assigning a gating value of 0.5 to a state in
a particular model may be used to simulate channel blocking by 50 percent.

%\smallskip \noindent \textbf{Rates} \\
Rates represent the kinetic rate constants which define the rate
of transition between two states. Rates are referred to by unique name. ModFossa
includes the following rate constant types: constant, Boltzmann voltage
dependence, exponential voltage dependence, and ligand-gated. The equations
for each rate constant are presented below.

\emph{Constant} rate constant are defined by a single parameter, $k$, like so:
\begin{equation}
 rate = k.
\end{equation}

\emph{Exponential} voltage dependent rate constants depend exponentially on
the membrane voltage: 
\begin{equation}
 rate(V) = a * exp(k*V),
\end{equation}
where $V$ is the membrane voltage, $exp$ is the exponential function, and $a$
and $k$ are parameters. 

\emph{Boltzmann} voltage dependence is defined using the sigmoid-like
Boltzmann equation:
\begin{equation}
 rate(V) = \frac{a}{\slfrac{1+exp[(V-V_{0.5})}{k}]},
\end{equation}
where $V$ is the membrane voltage, $V_{0.5}$ is the half-maximal activation
voltage (given as a parameter), $exp$ is the exponential function, and $a$
and $k$ are parameters. Taken from Angermann {\it et al.}\cite{angermann-2006}

\emph{Ligand-gated} rate constants depend on the concentration of a particular
ionic species, $S$, like so:
\begin{equation}
 rate([S]) = k * [S] ^ {n}, 
\end{equation}
where $[S]$ is the concentration of ligand $S$, $n$ is the ligand power, and
$k$ is a parameter.


%\smallskip \noindent \textbf{Connections} \\
Connections define a transition from one unique
state to another using a specified rate constant. States can have multiple
ingoing and outgoing connections. Rate constants may also be used in multiple
connections.
 
%\paragraph{Experiment definition:}
%\noindent
{\bf Experiment definition:}
An experiment consists of a channel Markov model, a voltage protocol, and an
optional concentration protocol. 
%\smallskip \noindent \textbf{Voltage protocol} \\
Voltage protocols are a fundamental technique in electrophysiology used to
measure the current response of ion channels. In the laboratory, a small
electrode is inserted into the cell, through the cell membrane. The electrode,
and therefore the cell membrane, is held constant at a desired voltage for a
length of time, as defined by the voltage protocol.

%\smallskip \noindent \textbf{Concentration protocol} \\
A concentration protocol defines the concentration values, extracellular
or intracellular, of a certain ligand. Like a voltage protocol, ModFossa's
concentration protocols are defined using \emph{hold stages}. However, there is
no support for stepped concentration protocol stages. During simulation, the
simulator will run the model through the entire voltage protocol using a single
stage from the concentration protocol. Next, the concentration protocol will
advance to the next stage, and the voltage protocol will be run using the new
concentration value. This continues until there are no stages left in the
concentration protocol.

%\paragraph{Data analysis and plotting:}   
{\bf Data analysis and plotting:}   
ModFossa supports the creation of a number of useful plots using the Python
plotting library, PyPlot. Unlike other simulators reviewed in Section 2, the
plots can be saved easily in a vector graphics format. Visual examples of each
plot type are found in Section 4.
%\smallskip \noindent \textbf{Model and experiment validation} \\
To assist the user in creating a valid simulation, ModFossa
provides validation methods which return detailed messages regarding any
problems with the model and experiment definitions.


\subsection{Implementation}
Sundials (SUite of Nonlinear and DIfferential/ALgebraic equation Solvers) is a C
library developed by Lawrence Livermore National Laboratory. ModFossa uses
Sundials to solve the system of ordinary differential equations governing the
ion channel state probabilities. Armadillo, a C++ linear algebra library
provides easy-to-use matrix structures and linear solvers, and is also used by
ModFossa.

\section{Results}
Several example models of varying complexity have been created and simulated
using ModFossa. A description of the each model is provided, along with a
listing of
the Python source code used define the model. Visual results are provided
in the form of plots generated by the code. Finally, an examination of
ModFossa's computational performance is provided.

\subsection{Simple model}
A simple three state model is presented in order to give the reader an
introduction to creating Markov models using ModFossa's Python interface.
This models has little biological relevance; for results
from a published calcium gated chloride ion channel model, see Section
\ref{section:angermann_model}. This section follows a tutorial format, with
explanations provided alongside a sequence of short source code snippets.


Let's create a simple model that will show some features we can
use to define real ion channel models, such as voltage steps, gating
variables, and voltage and ligand gated rate constants. Figure
\ref{fig:three_state_model} shows a Markov model with three states that we will
use for this example. $C_1$ is a closed state, while $O_1$ is a partially
conducting state with a gating variable of 0.5, and $O_2$ is a fully conducting,
open state. The rate constants for the model are provided in Table 2.
%\ref{table:three_state_model_rates}.

\begin{figure}
\begin{center}
%  \makebox[\textwidth]{\includegraphics[width=1.7in]{three_state_model.eps}}
\centerline{\includegraphics[width=1.7in]{three_state_model.eps}}
    \caption{Three state Markov model. The states 'C', 'O1', and 'O2' are
connected with the transition rates k1-k4.}
    \label{fig:three_state_model}
\end{center}
%\vspace*{-0.3in}
\end{figure}


%\begin{table}
%\begin{center}
% \tbl{Rate constants used in the three state model example.}
 % \label{table:three_state_model_rates}
 % \begin{tabular}{ l l l c }
  %  \hline
  %  Name & Type & Parameters & Units \\ \hline
  %  $k_1$ & ligand gated & $k=10e^{8}$, ligand = $Ca^{2+}$, power=1 &
%$M^{-1}s^{-1}$ \\
  %  $k_2$ & constant & k=1 & $s^{-1}$ \\
  %  $k_3$ & sigmoidal voltage gated & k=1, a=100, $V_{0.5}=50$ & $s^{-1}$ \\
  %  $k_4$ & constant & k=4 & $s^{-1}$ \\    
  %  \hline
 % \end{tabular}
%\end{center}
%\end{table}
%\begin{table}



Let's examine the code for the three state model. First, we define our states
and their gating variables. The closed state does not conduct, so we do not
have to specify its gating variable, it is 0 by default. State $O_1$ is a
partially conducting state, and $O_2$ is fully conducting.

{\small
\begin{singlespace}
\begin{lstlisting}
from modFossa import *
state('C')
state('O1', gating=0.5)
state('O2', conducting=True)
\end{lstlisting}
\end{singlespace}
}


Next, connect the states.

{\small
\begin{singlespace}
\begin{lstlisting}
connect(from_state='C', to_state='O1', rate='k1')
connect(from_state='O1', to_state='O2', rate='k2')
connect(from_state='O1', to_state='C', rate='k3')
connect(from_state='O2', to_state='O1', rate='k4')
\end{lstlisting}
\end{singlespace}
}

\begin{table}[h]
\vspace*{-0.3in}
\tbl{Rate constants used in the three state model example.}
  {\begin{tabular}{@{}cccc@{}} \toprule
    Name & Type & Parameters & Units \\ \colrule
    $k_1$ & ligand gated & $k=10e^{8}$, ligand = $Ca^{2+}$, power=1 &
$M^{-1}s^{-1}$ \\
    $k_2$ & constant & k=1 & $s^{-1}$ \\
    $k_3$ & sigmoidal voltage gated & k=1, a=100, $V_{0.5}=50$ & $s^{-1}$ \\
    $k_4$ & constant & k=4 & $s^{-1}$ \\    
    \botrule
  \end{tabular}}
\end{table}

Now we define our rate constants. We have a ligand gated rate constant which
depends
on the concentration of calcium. This could be any ligand that we chose, but we
must define its concentration in the experiment sweep, or ModFossa will give us
an error.
{\small
\begin{singlespace}
\begin{lstlisting}
rate('k1', type='ligandGated', k=10e8, ligand='Ca', power=1)
rate('k2', type='constant', k=1)
rate('k3', type='sigmoidal', k=1, a=100, v_half=50)
rate('k4', type='constant', k=4)
\end{lstlisting}
\end{singlespace}
}
Now, the remainder of the model definition.
{\small
\begin{singlespace}
\begin{lstlisting}
initialState('O2')
startAtSteadyState(False)
maxChannelConductance(1.16)
reversalPotential(0)
membraneCapacitance(100)
\end{lstlisting}
\end{singlespace}
}
Next, define the voltage protocol and experiment sweep. The voltage protocol
holds the voltage at -50 $mV$ for 500 $ms$, then depolarizes the membrane
by jumping 
to 50 $mV$ and holding for another half of a second. The concentration of
calcium is set to 500 $nM$ in the experiment sweep. This is necessary because
$k_1$ is a ligand gated rate dependent on $[Ca]$.

{\small
\begin{singlespace}
\begin{lstlisting}
voltageProtocol('vp')
voltageProtocolAddStage('vp', 'hold1', voltage=-50,  
            duration=500)
voltageProtocolAddStage('vp', 'hold2', voltage= 50,
            duration=500)
experimentSweep('three_state', 'vp', Ca=500e-9)
\end{lstlisting}
\end{singlespace}
}
Finally, validate and run the experiment and plot the results.

{\small
\begin{singlespace}
\begin{lstlisting}
validate()
run()
states = plotStates('three_state')
\end{lstlisting}
\end{singlespace}
}

Figure \ref{fig:three_state_probabilities} shows the results from the three
state model. Once again, we did not start in steady state, so we see the state
probabilities moving towards steady state under the
experiment sweep conditions of a -50 $mV$ membrane voltage, and 500 $nM$ $[Ca]$.
At 500 $mS$, the voltage protocol stage `hold2' is activated, and the voltage
is increased immediately to positive 100 $mV$. The effect of this jump in
voltage is evident in Figure \ref{fig:three_state_probabilities}. The
jump in the probability of state $O_1$ is due to a decrease in the value of the
sigmoidal rate constant, $k_3$.

%%%% Three state probabilities %%%%
\begin{figure}[h]
\begin{center} 
\makebox[\textwidth]{\includegraphics[width=3.5in]
{three_state_probabilities.eps}}
    \caption{A plot showing the evolution of the three state model state
probabilities as generated by ModFossa.}
    \label{fig:three_state_probabilities}
\end{center}
\vspace*{-0.3in}
\end{figure}

\subsection{Angermann CaCl currents model}
\label{section:angermann_model}

Our final example presents an implementation of the
calcium-activated chloride channel model developed by Angermann {\it et al.}\cite{angermann-2006} 
A brief description of the study is provided, along with
the model's description and implementation details. Later, the Python source
code used to implement the model in ModFossa is listed, along with several
figures and accompanying discussion. 

%\paragraph{Study details:}
{\bf Study details:}
In Ref.~\refcite{angermann-2006}, Angermann et al. attempt to determine how phosphatase
activity influences calcium-activated chloride channels in rabbit pulmonary
artery myocytes (smooth muscle cells). Whole cell patch clamp experiments were
conducted to measure the calcium-activated chloride current ($I_{ClCa}$)
under intracellular $[Ca^{2+}]$ ranging from 20 $mM$ to 1000 $mM$. The control
group contained 3 $mM$ ATP in the pipette solution, while the test group used 3
$mM$ adenosine $5'-(\beta, \gamma-imido)$-triphosphate (AMP-PNP).

Under both groups, the maximum $I_{ClCa}$ is greater at higher intracellular
$Ca^{2+}$ concentrations, as is expected for such channels. However, the
currents in test group (containing AMP-PNP) exhibited behavior indicative of a
negative shift in voltage dependent activation. The maximal conductance of
$I_{ClCa}$ was more than three-fold larger in the test group. Additionally, the
activation times were less and the deactivation time greater in the test group
compared to the control. These observations led the authors to the conclusion
that the increased $I_{ClCa}$ activity in the test group (pipette
containing AMP-PNP) is not explained by $Ca^{2+}$ channel sensitivity, nor the
number of $Ca^{2+}$ activating the channel, but instead by a negative shift in
the voltage dependent activation. 

%\paragraph{Model description:}
{\bf Model description:}
To simulate the behavior of the test and control groups, Angermann et al.
created a Markov chain kinetic model based on the work by Kuruma and Hartzell.\cite{kuruma-2000} 
Slightly different parameters were chosen for the test and
control groups in an attempt to make clear the difference in behavior between
the two. The model, shown in Figure \ref{fig:angermann_model}, consists of four
closed states and three open states. Calcium binding occurs in the closed states
with rates directly proportional to $[Ca]$. State $C_1$ is representative of the
channel's state when no $Ca^{2+}$ ions are bound to the receptor sites, while
state $C_4$ represents the channel with all three $Ca^{2+}$ binding sites
occupied. The closed states with at least one binding site occupied can
transition into open states with constant channel opening rates. Closed states
with more occupied binding sites have higher values for the channel opening
rates, giving the model its concentration-dependent behavior. Voltage
dependent channel closing rates take the model from open to closed states. 

In Ref.~\refcite{angermann-2006}, the behavior of the channel in the presence of
AMP-PNP was simulated by setting the gating variables of all three of the
channel's open states to 1. Conversely, to simulate the channel in the presence
of ATP, only channel $O_1$ was assigned a gating variable of 1, while the rest
of the open states were given a value of zero. This reproduces the
blocking effect by phosphorylation that the authors hypothesize as a possible
explanation for the inhibition of $I_{ClCa}$ while the channel is subjected to
ATP. In addition to reducing the channel gating variables to model the effect of
ATP, the authors also increased the magnitude of the channel closing rates and
shifted the voltage dependence towards more positive potentials.

By altering the Markov model's parameters in an insightful manner, Angermann et
al. were able to simulate the behavior of $Ca^{2+}$-activated $Cl^-$ currents
in the presence of both AMP-PNP and ATP. 

%%%% Markov Model %%%%
\begin{figure}[]
\begin{center}
\makebox[\textwidth]{\includegraphics[width=3in]
{angermann_model.eps}}
    \caption{Angermann CaCl Markov model. The model has four closed states and
three open states. State C2 represents the channel with a single $Ca^{2+}$
binding site occupied, while state C4 represents the channel with all three
binding sites occupied. Note that the model can only transition to an open
state in the presence of at least one bound $Ca^{2+}$.}
    \label{fig:angermann_model}
\end{center}
\end{figure}

%\paragraph{Model implementation in ModFossa:}
{\bf Model implementation in ModFossa:}
The calcium activated chloride channel Markov model shown in Figure
\ref{fig:angermann_model} was implemented in ModFossa using parameter values 
chosen to match the behavior of $I_{ClCa}$ in the
presence of ATP, and were modified slightly from Thibeault
{\it et al.}\cite{angermann-2006}

%\paragraph{ModFossa output:}
{\bf ModFossa output:}
ModFossa can generate a number of IV curves at a given time, and place
them on a single plot, as shown in Figure \ref{fig:angermann_ca_late}. This
figure shows the \emph{late IV} curves, measured at 1099 milliseconds into the
simulation, which is the moment of maximum channel current. 

To show the effects of both ligand concentration and voltage on the channel
conductance, ModFossa's \emph{plotGvsConcentration} method can be used to
generate a chord conductance plot (not shown). This plot contains 13 traces of
conductance as a function of the internal calcium concentration. Each trace
corresponds to a different step value in the voltage protocol.

A second chord conductance plot was generated using ModFossa (Figure
\ref{fig:angermann_ca_g_vs_v}). This plot shows channel conductance as a
function of voltage, for 6 different concentration values. The simulated data
(Figure \ref{fig:angermann_ca_g_vs_v}) match the experimental data well.


%%%% Late IV %%%%
\begin{figure}[]
\begin{center}
  \makebox[\textwidth]{\includegraphics[width=5in]{angermann_ca_late.eps}}
    \caption{Angermann CaCl late IV curves given various internal $Ca^{+2}$
concentrations. This figure was generated using ModFossa's \emph{plotIV}
method at time 1099 milliseconds.}
    \label{fig:angermann_ca_late}
\end{center}
\end{figure}

%%%% G vs V %%%%
\begin{figure}[]
\begin{center}
 
\makebox[\textwidth]{\includegraphics[width=4in]
{angermann_ca_g_vs_v.eps}}
    \caption{Simulated Angermann CaCl conductance at varying concentrations, as
a function of membrane voltage. This figure was
generated using ModFossa's \emph{plotGvsV} method.}
    \label{fig:angermann_ca_g_vs_v}
\end{center}
\end{figure}


\section{Discussion}

\subsection{Summary}
We have identified and addressed a gap in ion channel simulation software by
developing ModFossa, a Python library dedicated to the construction and
simulation of ion channels based on continuous time Markov processes. While
several recent software tools have been released for this purpose, none are
Python libraries. Python is a popular tool among research scientists, therefore
we hypothesize that a versatile and easy-to-use ion channel simulation library
has the potential for providing a valuable service to those developing and
testing Markov models of ion channels. In summary, ModFossa allows for the
creation of channel models using simple Python syntax. Several rate constant
types are supported, including exponential and Boltzmann voltage dependent
rates, as well as ligand-gated rates. Voltage and current protocols are defined
easily. Some common plots including IV curves, state probabilities, conductance
versus voltage, and conductance versus ligand concentration are generated using
a single function call. Finally, ModFossa’s core is implemented in C++ using the
Sundials  solver, providing the benefit of fast execution
times. 

\subsection{Applications}
ModFossa can provide rapid model development and testing for researchers with
varying programming skill. A common task in ion channel model development is
parameter searching, which requires the modeler to find the rate constant
parameters which best fit the experimental data. Because ModFossa is usable from
Python, existing machine learning libraries such as scikit-learn can be
leveraged to perform parameter searches automatically. This situation is where
ModFossa’s fast execution speed provides a great benefit, as it is not uncommon
for modelers to leave parameter searches running for several days, or weeks. 

\subsection{Future work}
Several ideas for additional development of ModFossa are presented here.
First, some minor enhancements such as the addition of more sample models,
improved customization of plotting functions (plot size, legend, colors, data
range, etc.), and standardization of the data structures used in the results
module would add value to the software by ensuring a smooth experience for new
users. Additionally, creating and distributing a Debian package for ModFossa
would ease installation significantly. 

Several preliminary users have expressed interest in user-defined rate contant
equations. A custom rate equation could be written by the user in Python using
either a decorator function, or by inheriting from a base class. During
simulation, the simulator would call the user's custom rate equation, passing in
the current state of the simulation (membrane voltagte, concentrations, time,
etc.), and use the returned value as the rate for whichever transitions the user
specified. The performance hit due to calling a Python function from C++ during
the simulation would have to be quantified, but it is not expected to be
significant.

Lastly, the idea of offering ModFossa as a web service with a rich graphical
user interface has been explored. Doing so would encourage use from a wider
audience, particularly from those with limited or no programming experience.

%\clearpage

\section*{Acknowledgments}
This material is based in part upon work supported by: The National Science Foundation under grant number(s) IIA-1329469. 
Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not 
necessarily reflect the views of the National Science Foundation.

\bibliographystyle{ws-jbcb}
\bibliography{bibliography-1}


%\appendix
% \vspace*{13pt}
 \noindent
 \parbox{5truein}{
 \begin{minipage}[b]{1truein}
 \centerline{{\psfig{file=Garreth.jpg,width=1in}}}
 \end{minipage} 
 \hfill         %to 2nd column
 \begin{minipage}[b]{3.85truein}
 {\bf Garreth B. Ferneyhough} received a B.S. and an M.S. in 
Computer Science and Engineering from the University of Nevada, Reno
in 2010 and 2013 respectively.  His research is in the area of Computational
Neuroscience. He is currently working as a Software Engineer in the private sector.
\vspace*{0.45in}
 \end{minipage} } 




%-------------------------------------------------
 \vspace*{13pt}	%Corey Thibeault
 \noindent
 \parbox{5truein}{
 \begin{minipage}[t]{1truein}
 \centerline{{\psfig{file=TM-Corey.jpg,width=1in}}}
 \end{minipage} 
 \hfill         %to 2nd column
 \begin{minipage}[b]{3.85truein}
 {{\bf Corey M. Thibeault} received a B.S.
degree in Computer Engineering from
the Rochester Institute of Technology,
Rochester, NY, USA, in 2004, the
B.S. and M.S. degrees in Mechanical
Engineering from Minnesota State
University, Mankato, MN, USA, in 2006
and 2009, respectively, and the M.S.
degree in Computer Engineering and the
Ph.D. degree in Biomedical Engineering
from the University of Nevada, Reno, NV, 
in 2012. He is currently a Senior Data Scientist 
at\hfilneg} 
 \end{minipage} 

%\vspace*{4.8pt} 
 \noindent 
 
Neural Analytics, Los Angeles, CA, USA
}

%\eject

%-------------------------------------------------

 \vspace*{13pt}	%Sergiu Dascalu
 \noindent
 \parbox{5truein}{
 \begin{minipage}[t]{1truein}
% \centerline{{\psfig{file=dascalu.jpg,width=1in,height=1.28in}}}
\centerline{{\psfig{file=dascalu-2.jpg,width=1in}}} \end{minipage} 
 \hfill         %to 2nd column
 \begin{minipage}[b]{3.85truein}
 {{\bf Sergiu M. Dascalu} is a Professor in the Department of Computer
Science and Engineering at the University of Nevada,
Reno, USA, which he joined in 2002. In 1982 he received a
Master’s degree in Automatic Control and Computers from
the Polytechnic University of Bucharest, Romania and in
2001 a Ph.D. in Computer Science from Dalhousie University,
Halifax, NS, Canada. His main research interests are
in the areas of software engineering and \hfilneg} 
 \end{minipage} 

%\vspace*{4.8pt} 
 \noindent
human–computer
interaction. He has published  over 140 peer-reviewed papers
and has been involved in numerous projects funded
by industrial companies as well as federal agencies such as
NSF, NASA, and ONR.}

%\eject
%---------------------------------------------------------

 \vspace*{13pt}	%Fred Harris
 \noindent
 \parbox{5truein}{
 \begin{minipage}[t]{1truein}
 \centerline{{\psfig{file=harris.jpg,width=1in,height=1.28in}}}
 \end{minipage} 
 \hfill         %to 2nd column
 \begin{minipage}[b]{3.85truein}
 {{\bf Frederick C. Harris, Jr.}. is a Professor in the Department
of Computer Science and Engineering and the Director of
the High Performance Computation and Visualization Lab
and the Brain Computation Lab at the University of Nevada,
Reno, USA. He received his B.S. and M.S. degrees in Mathematics
and Educational Administration from Bob Jones
University in 1986 and 1988 respectively, and his M.S. and
Ph.D. degrees in Computer Science from Clemson University
in 1991 and 1994 respectively.\hfilneg} 
 \end{minipage} 

%\vspace*{4.8pt} 
 \noindent
 He is a SeniorMember of
ACM and ISCA. His research interests
are in parallel computation, computational neuroscience,
computer graphics and virtual reality.}

 \vfill\eject

\end{document}

