\graphicspath{ {./images/} }

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.

\section{Simple models}
Simple two and three state models are presented in order to give the reader an
introduction to creating Markov models using ModFossa's Python interface.
The simple models in this section have 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.


\subsection{Two State Model}
The first example provided is a simple ion channel model with two states, Open,
and Closed (Figure \ref{fig:two_state_model}). The model defines a transition
from
the Open state to the Closed state with a rate $k_{off}$ of 1 $s^{-1}$, and a
transition from the Closed state to the Open state with a rate $k_{on}$ of
10 $s^{-1}$. We define the initial state of the model as the Closed state.

%%%% Two state model %%%%
\begin{figure}[h!]
\begin{center}
  \makebox[\textwidth]{\includegraphics[width=3in]{two_state_model.eps}}
    \caption{A two state Markov model containing two states, 'Open', and
'Closed.' The arrows between the states represent transition rates.}
    \label{fig:two_state_model}
\end{center}
\end{figure}

The two state model is implemented using ModFossa's Python interface. Let's
examine the source code in detail. The first step is to import the library:

\begin{singlespace}
\begin{python}
from modFossa import *
\end{python}
\end{singlespace}


Then we define our two states and name them `Open' and `Closed', and make
`Open' a conducting state. Note that although we are defining `Open' as
conducting, it will have no effect on the results presented for the two state
model, because we are not calculating the channel current or conductance, only
state probabilities.

\begin{singlespace}
\begin{python}
state('Open', conducting=True)
state('Closed')
\end{python}
\end{singlespace}

Next, the connections between the two states are created. ModFossa allows
the user to define connections before the states or rate constants have
been defined.

\begin{singlespace}
\begin{python}
connect(from_state='Closed', to_state='Open', rate='kon')
connect(from_state='Open', to_state='Closed', rate='koff')
\end{python}
\end{singlespace}


Now we define the rate constants, $k_{on}$ and $k_{off}$, both of type
`constant'. 
\begin{singlespace}
\begin{python}
rate('kon', type='constant', k=10)
rate('koff', type='constant', k=1)
\end{python}
\end{singlespace}

The last piece of information needed to describe the behavior of the model is
the initial state, so we set it to closed. We also instruct ModFossa not
to set the initial conditions of the simulation to the steady state
probabilities, and to instead use the initial state, `Closed' as the initial
conditions. This allows us to observe the state probabilities of `Open' and
`Closed' approach steady state, as shown in Figure
\ref{fig:two_state_probabilities}.

\begin{singlespace}
\begin{python}
initialState('Closed')
startAtSteadyState(False)
\end{python}
\end{singlespace}

Although the simple two state model has no use for them, we must set the
maximum channel conductance, and the reversal potential. These are required for
the simulation to run. 

\begin{singlespace}
\begin{python}
maxChannelConductance(0) #mS/cm^2
reversalPotential(0)     #mV
\end{python}
\end{singlespace}

Now we define our voltage protocol. We give it a name, `vp', and add a single
stage to it with the name `hold'. The `hold' stage will hold the membrane
potential at -50 $mV$ for 500 $ms$. Next, an experiment sweep with the name
`twoStateModelSweep' is created. It used the voltage protocol we just defined,
`vp'. In later examples we will use the experiment sweep to define the
ionic concentrations in addition to the voltage protocol. 

\begin{singlespace}
\begin{python}
voltageProtocol('vp')
voltageProtocolAddStage('vp', 'hold', 
  voltage=-50, duration=500)
experimentSweep('twoStateModelSweep', 'vp')
\end{python}
\end{singlespace}

Finally, we can validate our model, run our experiment sweeps, and plot the
results. We have only defined one experiment sweep, but if we had more they
would run as well.

\begin{singlespace}
\begin{python}
validate()
run()
plotStates('twoStateModelSweep')
\end{python}
\end{singlespace}

Figure \ref{fig:two_state_probabilities} shows the evolution of the state
probabilities as a plot generated by ModFossa. At time
$t=0 ms$, the model is the Closed state. As
time advances, the probabilities move towards steady state. At $t=500 ms$, the
probability of the Open state (rounded to three decimal places) is 0.905.
Let's compare that to the analytical solution of the steady-state probability.

The differential equation describing the Open state probability is

\begin{equation}
\frac{dP_{Open}}{dt} = k_{on}(1-P_{Open})-k_{off}P_{Open}. 
\end{equation}

Steady-state occurs when the state probabilities are not changing, i.e.,
$\frac{dP_{Open}}{dt} = 0$. The steady-state Open probability is therefore

\begin{equation}
P_{OpenSS} = \frac{k_{on}}{k_{on}+k_{off}}.
\end{equation}

Substituting in our rate constants from the two state model gives us
\begin{equation}
P_{OpenSS} = \frac{10}{10 + 1} = 0.909.
\end{equation}

So, after 500 $ms$ of simulation time, we are within $99.6\%$ of the true
steady-state probabilities.

%%%% Two state probabilities %%%%
\begin{figure}[h!]
\begin{center}
  \makebox[\textwidth]{\includegraphics[width=6in]{two_state_probabilities.eps}}
    \caption{Evolution of the `Open' and `Closed' state probabilities of the
simple two state model. At the initial conditions at time t=0, the system is
entirely in the Closed state. As the simulation runs, the probabilities
approach steady state.}
    \label{fig:two_state_probabilities}
\end{center}
\end{figure}

\subsection{Three State Model}
Now let's create a more complex model that will show some features we can
actually 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
\ref{table:three_state_model_rates}.

\begin{figure}[h!]
\begin{center}
  \makebox[\textwidth]{\includegraphics[width=3in]{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}
\end{figure}

\begin{table}[h!]
\begin{center}
  \caption{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}

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.

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

Next, connect the states.

\begin{singlespace}
\begin{python}
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{python}
\end{singlespace}

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.
\begin{singlespace}
\begin{python}
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{python}
\end{singlespace}

We must supply the remainder of the model definition. This is similar to the
previous two state model example.

\begin{singlespace}
\begin{python}
initialState('O2')
startAtSteadyState(False)
maxChannelConductance(1.16)
reversalPotential(0)
membraneCapacitance(100)
\end{python}
\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]$.

\begin{singlespace}
\begin{python}
voltageProtocol('vp')
voltageProtocolAddStage('vp', 'hold1', voltage=-50,  
			duration=500)
voltageProtocolAddStage('vp', 'hold2', voltage= 50,
			duration=500)
experimentSweep('three_state', 'vp', Ca=500e-9)
\end{python}
\end{singlespace}

Finally, validate and run the experiment and plot the results.


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

Figure \ref{fig:three_state_probabilities} shows the results from the three
state model. One 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=6in]
{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}
\end{figure}

%\clearpage
\section{Angermann CaCl currents model}
\label{section:angermann_model}

Our final example presents an implementation of the
calcium-activated chloride channel model developed by Angermann 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. 

\subsection{Study details}
In \cite{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. 

\subsection{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 \cite{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}[h!]
\begin{center}
\makebox[\textwidth]{\includegraphics[width=5in]
{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}

\subsection{Model implementation in ModFossa}

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

\subsection{Code listing}
The Python source code for the model is provided on the following page. 
\begin{singlespace}
\pythonexternal{angermann_cacl_code.py}
\end{singlespace}

\cleardoublepage

%%%% Model Parameters Table %%%%
\begin{table}
\begin{center}
  \caption{Parameters for Angermann CaCl model in presence of ATP. Reproduced
with new values from \cite{thibeault-2011}.}
  \label{table:angermann_parameters}
  \begin{tabular}{ l c c }
    \hline
      & Value & Units \\ \hline
    Conductance and & & \\
    equilibrium potential & & \\
    Maximal conductance & 1.16 & $\slfrac{mS}{cm^2}$ \\ 
    $E_{Cl}$ & 0 & $mV$ \\
    Values of gating variables & & \\
    $C_1$ & 0 & \\ 
    $C_2.Ca$ & 0 & \\ 
    $C_3.2Ca$ & 0 & \\ 
    $C_4.3Ca$ & 0 & \\ 
    $O_1$ & 1 & \\ 
    $O_2$ & 0.5 & \\ 
    $O_3$ & 0 & \\  
    $Ca^{2+}$ binding rates & & \\
    $k_{on}(C_1 \rightarrow C_2.Ca)$ & $20 \times 10^6$ & $M^{-1}s^{-1}$ \\
    $k_{on}(C_2.Ca \rightarrow C_3.2Ca)$ & $20 \times 10^6$ & $M^{-1}s^{-1}$ \\
    $k_{on}(C_3.2Ca \rightarrow C_4.3Ca)$ & $20 \times 10^6$ & $M^{-1}s^{-1}$
\\ 
    Unbinding rates & & \\
    $k_{off}(C_2.Ca \rightarrow C_1 + Ca)$ & $50$ & $s^{-1}$ \\
    $k_{off}(C_3.2Ca \rightarrow C_2.Ca + Ca)$ & $50$ & $s^{-1}$ \\
    $k_{off}(C_4.3Ca \rightarrow C_3.2Ca + Ca)$ & $50$ & $s^{-1}$ \\
    Channel opening rates & & \\
    $\alpha_{1}(C_2.Ca \rightarrow O_1)$ & $1$ & $s^{-1}$ \\    
    $\alpha_{2}(C_3.2Ca \rightarrow O_2)$ & $25$ & $s^{-1}$ \\ 
    $\alpha_{3}(C_4.3Ca \rightarrow O_3)$ & $200$ & $s^{-1}$ \\  
    Channel closing rates & & \\
    $\beta_1(V)(O_1 \rightarrow C_2.Ca)$ & & \\
    where $a, V_{0.5},$ and $k =$ & $ 60, -40,$ and $40$ & $s^{-1}, mV, mV$
\\
    $\beta_2(V)(O_2 \rightarrow C_3.2Ca)$ & & \\
    where $a, V_{0.5},$ and $k = $ & $ 35, 0,$ and $50$ & $s^{-1}, mV, mV$
\\ 
    $\beta_3(V)(O_3 \rightarrow C_4.3Ca)$ & & \\
    where $a, V_{0.5},$ and $k = $ & $ 25, 140,$ and $40$ & $s^{-1}, mV, mV$
\\ \hline
    \multicolumn{3}{l}{All $\beta_x(V)$ use the following Boltzmann equation
form:} \\
    \multicolumn{3}{l}{ $\beta_x(V) =
\slfrac{a}{\slfrac{\{1+exp[(V-V_{0.5})}{k}]\}}$.} \\
	    
  \end{tabular}
\end{center}
\end{table}
\clearpage




  
\subsection{ModFossa output}
Figure \ref{fig:angermann_ca_currents} shows the figure generated by ModFossa
using the \emph{plotCurrents} method. Each current trace shows the
current passing through the channel at various $Ca^{+2}$ levels. As expected,
higher $[Ca^{+2}]$ invoke larger channel currents. The bottom item in the
figure is the voltage protocol used in the experiment. Note that labels and
calibration bars were added manually using Inkscape. 

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. Similarly, Figure
\ref{fig:angermann_ca_tail} shows the \emph{tail IV} curves, which were
measured at 1100 milliseconds.  

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, as shown in
Figure \ref{fig:angermann_ca_g_vs_cai}. 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. As shown in
Figure \ref{fig:angermann_ca_currents}, the voltage protocol contains 13
stepped stages, hence the 13 traces in Figure \ref{fig:angermann_ca_g_vs_cai}.

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. A chord
conductance plot of Angermann et. al's measured experimental data is shown in
Figure \ref{fig:angermann_g_vs_v_atp_experimental}. As shown, the simulated data
(Figure \ref{fig:angermann_ca_g_vs_v}) match the experimental data (Figure
\ref{fig:angermann_g_vs_v_atp_experimental}) well.

%%%% Currents %%%%
\begin{figure}[h!]
\begin{center}
 
\makebox[\textwidth]{\includegraphics[width=3.2in]
{angermann_ca_currents_single_column_labeled.eps} }
    \caption{Figure generated by ModFossa showing Angermann CaCl currents given
various internal $Ca^{+2}$ concentrations. The voltage protocol is shown at the
bottom of the figure.}
    \label{fig:angermann_ca_currents}
\end{center}
\end{figure}

%%%% Late IV %%%%
\begin{figure}[h!]
\begin{center}
  \makebox[\textwidth]{\includegraphics[width=7in]{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}

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

%%%% G vs Cai %%%%
\begin{figure}[h!]
\begin{center}
\makebox[\textwidth]{\includegraphics[width=5in]
{angermann_ca_g_vs_cai.eps} }
    \caption{Angermann CaCl conductance at varying voltages as a function of
intracellular $Ca^{+2}$ concentration. There are 13 traces on the plot, which
correspond to the 13 different steps in the voltage protocol. This figure was
generated using ModFossa's \emph{plotGvsConcentration} method.}
    \label{fig:angermann_ca_g_vs_cai}
\end{center}
\end{figure}

%%%% G vs V %%%%
\begin{figure}[h!]
\begin{center}
 
\makebox[\textwidth]{\includegraphics[width=5in]
{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}

%%%% G vs V Experimental%%%%
\begin{figure}[h!]
\begin{center}
 
\makebox[\textwidth]{\includegraphics[width=5in]
{angermann_g_vs_v_atp_experimental.eps}}
    \caption{Measured data from \cite{angermann-2006} showing conductance as a
function of voltage and concentrations. The data were fit to a Boltzmann
function using the least squares method.}
    \label{fig:angermann_g_vs_v_atp_experimental}
\end{center}
\end{figure}

\subsection{Runtime Analysis}
The runtimes of ModFossa versus a similar Matlab implementation are provided in
Figure \ref{fig:runtime_graph}. ModFossa is approximately 17 times faster.

%%%% Runtime %%%%
\begin{figure}[h!]
\begin{center}
 
\makebox[\textwidth]{\includegraphics[width=4in]{runtime.eps}}
    \caption{ModFossa runtime versus a Matlab implementation. The model and
voltage protocol used are described in Section \ref{section:angermann_model}.
The x-axis corresponds to the number of experiment sweeps performed at each data
point.}
    \label{fig:runtime_graph}
\end{center}
\end{figure}


