

This section discusses two biomathematical models
that have been analyzed using techniques
outside of the machine learning field.
The models will be introduced, along with a description
of the experimental goals and results of the experiments
using the SPY framework.

\subsection{Models Under Consideration}\label{sect:biomath}

\subsubsection{Competitive Exclusion}
\cite{biomodeling80} develops dozens of biomathematical
models.  
Dynamic population models are one of the classic types of 
biomathematical models.
Populations are modeled using coupled differential
equations.  The first derivatives (in time) are  
constructed from first principles.
These derivatives usually contain several parameters.
%that,
%while they have some physical meaning, are too abstract
%to be measured by field biologists
These parameters can be studied analytically to
define different types of behaviors of the model.
We will be studying a system of 2 different species
that compete in some way in the same niche.
We can think of this competition in terms of
food supply, space, toxicity, or anything
else that would lower the carrying capacity
of the niche for both species.
This system is called competitive exclusion
because analytic analysis shows that under most
conditions, one species will be driven to
extinction.
The general form of the population equations is
\begin{align}
\frac{dN_1}{dt} = r_1 N_1 
	\left( 1 - \frac{N_1}{K_1} - b_{12} \frac{N_2}{K_1} \right)
        \label{eqn:n1}\\
\frac{dN_2}{dt} = r_2 N_2 
	\left( 1 - \frac{N_2}{K_2} - b_{21} \frac{N_1}{K_2} \right)
        \label{eqn:n2}
\end{align}
where $N_1$ and $N_2$ are the size of the two species,
$K_1$ and $K_2$ are the carrying capacities for
the two species in the absence of the other species,
$r_1$ and $r_2$ are the growth rates of the 
species, and $b_{12}$ and $b_{21}$ measure
the degree to which the two species effect each other.
All these constants are positive and real-valued.
Following the methodology of \cite{biomodeling80}
we rewrite Equations \ref{eqn:n1} and \ref{eqn:n2} 
in dimensionless terms as
\begin{align}
\frac{du_1}{d\tau} = u_1(1 - u_1 - a_{12}u_2) \label{eqn:compEx1} \\
\frac{du_2}{d\tau} = \rho u_2(1 - u_2 - a_{21}u_1) \label{eqn:compEx2}
\end{align}
where $u_1=\frac{N_1}{K_1}$, $u_2=\frac{N_2}{K_2}$,
$a_{12}=b_{12}\frac{K_2}{K_1}$, 
$a_{21}=b_{21}\frac{K_1}{K_2}$, $\rho = \frac{r_2}{r_1}$,
and $\tau=r_1t$.
We have eliminated one term (we scaled $r_1$ to one
and redefined $r_2$ to a ratio) and normalized
our population sizes so that the carrying capacity
of each species is 1.


\subsubsection{Animal Neurons}
A simple model of animal neurons was
developed in \cite{izhikevich03}.
This model has the advantages of being
computationally efficient, while also
being able to simulate most of the
behaviors of much more complex models.
Its main drawback is that its parameters have
no physical significance.
The model is a pair of coupled
ordinary differential equations.
The entire model has the form
\begin{align}
\frac{dv}{dt} &= 0.04v^2 + 5v + 140 - u + I \label{eqn:neuron1}\\
\frac{du}{dt} &= a(bv - u) \\
&\text{if} ~v \geq 30 mV, \text{then} ~v \leftarrow c~\text{and}~u \leftarrow u + d \notag
\end{align}
Where $v$ and $u$ are dimensionless variables,
$a,b,c,d$ are dimensionless parameters,
and $I$ is an input DC signal.
The variable $v$ represents the membrane potential
of the neuron and is the experimentally observable quantity
in this model.
The variable $u$ is a recovery variable, which models
the rate at which $K^+$ and $Na^+$ channels can open and close
in the neuron.

\cite{izhikevich03} and \cite{izhikevich04} state that
almost all behaviors seen in real neurons can be modeled
using this model, depending
on the settings of the different parameters.  
There is matlab program available
for download\footnote{
\url{http://vesicle.nsi.edu/users/izhikevich/publications/izhikevich.m}
} that lists the different values
of the parameters that lead to different
classes of behavior.

\subsection{Previous Analytic Work with Models}\label{sect:analysis}
\subsubsection{Competitive Exclusion} \label{sect:compExAnalysis}
The dynamic population model we are using
is simple enough to be analyzed using a number
of different techniques.
Here we simply summarize the 
linear algebra based technique used
in \cite{biomodeling80}.
First we construct the matrix $A$ around
a point where 
$\frac{du_1}{d\tau} = \frac{du_2}{d\tau} = 0$.
Then the sign of the eigenvalues
of the matrix $A$ determine whether
the equilibrium points are stable or unstable.
After finding the eigenvalues at the four
different equilibrium points, it turns out
there is only a single set of conditions
that will lead to a stable equilibrium
point where $N_1,N_2 \neq 0$.
If $0 \leq a_{12} \leq 1$ and $0 \leq a_{21} \leq 1$ then there will
be a stable equilibrium point at
$(u^*_1, u^*_2)$, where 
$0 < u^*_1 < 1$ and $0 < u^*_2 < 1$ 
(recall equations \ref{eqn:compEx1} and \ref{eqn:compEx2}
were nondimensionalized, so the carrying capacity
for the two different species is 1).


\subsubsection{Animal Neurons}
Our neuron model is much harder to analyze
symbolically, both because it is not first-order
(recall the $v^2$ term in equation \ref{eqn:neuron1})
and because of the reset condition.
Izhikevich has spent much time describing
some of the analytic properties of different models 
of animal neurons \cite{izhikevich2000}.
All of these discussion are beyond the scope of this thesis.  
It should be noted that this
analysis takes dozens of pages and a 
high degree of mathematical training to discover.
The various parameter settings that lead to
the different behaviors are published,
but no explanation is given as to how
they were found (which suggests a 
simple trial-and-error search).

\subsection{Experimental Goals}
There are two slightly different goals we have
in mind for our experiments.
We would like our framework to find the relevant
range
restrictions and we would like it to NOT offer
any range restrictions on parameters that
are irrelevant.

\subsubsection{Competitive Exclusion}
For the dynamic population model we
want to discover the conditions that 
create a stable equilibrium point where
$N_1,N_2 \neq 0$, as discussed in 
\ref{sect:compExAnalysis}.
The range restrictions that lead
to this behavior are $0 \leq a_{12} \leq 1$ and $0 \leq a_{21} \leq 1$.
Since symbolic analysis shows that the value of 
$\rho$ and the initial populations do not affect
the location or existence of equilibrium points,
we hope that SPY will not offer any range restrictions
on those variables.


\subsubsection{Animal Neurons}
Of all the different behaviors described in 
\cite{izhikevich04}, three were chosen
because it was easy to
describe them in mathematical terms.
First, {\it phasic spiking}, \fig{ps}, which is a single spike and reset.
Second, {\it tonic bursting}, \fig{tb}, 
which is several spikes close together in time, 
followed by a quiescent state, repeated several times.
Third, {\it phasic bursting}, \fig{pb}, 
is several spikes close together in time, 
followed by a quiescent state without any repeats.
We attempted to find the input value ranges
that lead to these different behaviors.
\begin{figure}
\centering
\subfloat[phasic spiking]{\label{fig:ps}{\includegraphics[width=1.71in]{Figures/phasic-spiking}}}
\subfloat[tonic bursting]{\label{fig:tb}{\includegraphics[width=1.71in]{Figures/tonic-bursting}}}
\subfloat[phasic bursting]{\label{fig:pb}{\includegraphics[width=1.71in]{Figures/phasic-bursting}}}
\caption{Behavior modes of animal neuron}\label{fig:phasicSpiking}
\end{figure}


\subsection{Methodology}

\subsubsection{Competitive exclusion model in the SPY language}
To highlight some of the interesting features of the
SPY language we will discuss the particulars
of the competitive exclusion model, as it was
developed in SPY.
The SPY code for the competitive exclusion model is
listed in \fig{varPicking} - \fig{execution}.

\begin{figure}
\begin{footnotesize}
\begin{listing}[5]{1}
Choice(a1SpyVal,=,n){
  return linear(a1SpyVal_min,a1SpyVal_max,n); }
Choice(k1SpyVal,=,n) {
  return linear(k1SpyVal_min, k1SpyVal_max,n); }
Choice(k2SpyVal,=,n) {
  return linear(k2SpyVal_min, k2SpyVal_max,n); }
Choice(alphaSpyVal,=,n){
  return linear(alphaSpyVal_min,alphaSpyVal_max,n); }
Choice(betaSpyVal,=,n) {
  return linear(betaSpyVal_min,betaSpyVal_max,n); }

function rho(   tmp) {
  tmp = a1SpyVal();
  #split 0-1 in half and then rescale
  if(tmp>.5) return 28*tmp-13;
  else return tmp*2;
}
function a12(   tmp){
  tmp = alphaSpyVal();
  #split 0-1 in half and then rescale to 0-1 and 1-15
  if(tmp>.5) return 28*tmp-13;
  else return tmp*2;
}
function a21(   tmp) {
  tmp = betaSpyVal();
  #split 0-1 in half and then rescale to 0-1 and 1-15
  if(tmp>.5) return 28*tmp-13;
  else return tmp*2;
}
# don't allow initial populations of zero
function X0() {
  return (k1SpyVal()+.02)*2;
}
function Y0() {
  return (k2SpyVal()+.02)*2;
} 
\end{listing}
\end{footnotesize}
\caption{Picking model variables in the competitive exclusion model}\label{fig:varPicking}
\end{figure}


\begin{figure}
\begin{footnotesize}
\begin{listing}[4]{1}
# Model:
#    dx/dt = x * (1 - x - a12*y)
#    dy/dt = rho*y * (1 - y - a21*x)

function dx(x,y) {
  return x*(1-x-a12()*y);
}
function dy(x,y) {
  return rho()*y*(1-y-a21()*x);
}
\end{listing}
\end{footnotesize}
\caption{Competitive exclusion model, adapted from\cite{biomodeling80}}
\label{fig:model}
\end{figure}

\begin{figure}
\begin{footnotesize}
\begin{listing}[4]{1}
function main(warmup,   n,delx,dely,xt,yt) {
  a1SpyVal(n);                     
  k1SpyVal(n);k2SpyVal(n);
  alphaSpyVal(n);betaSpyVal(n);    

  xt=X0();yt=Y0();
  REDO=0
  if (!warmup) {print "# new experiment" >> "points.dat";
    print count++ " " a12() " " a21() >> "parameters.dat";}
  do {
    REDO++;
    if (!warmup) print xt " " yt " " sumworth >> "points.dat";

    # calculate change to populations
    delx = dx(xt,yt)/UNITSTEP;
    dely = dy(xt,yt)/UNITSTEP;
    xt += delx;xt=(xt<SMALL?0:xt);
    yt += dely;yt=(yt<SMALL?0:yt);
    totalChange = (abs(delx)<SMALL?0:delx) + (abs(dely)<SMALL?0:dely);
  } while (totalChange != 0 || REDO<10)
  return ((xt>.SMALL?1:0)+(yt>.SMALL?1:0));
}
\end{listing}
\end{footnotesize}
\caption{Model execution code in the competitive exclusion model}\label{fig:execution}
\end{figure}


\fig{varPicking} shows the use of the {\tt Choice} function.
This function is the main interface between the model
and the learner.
The {\tt Choice} function maintains a cache of
all the values it has returned during an iteration,
and these values are what the learner takes as training data.
At the start of each run through the model, the function
picks a new value and always returns the same value
until the start of another run.
Currently SPY interacts with its learner best if the values
remembered by the {\tt Choice} function are all
$[0-1]$, hence the five functions after the 
{\tt Choice} functions in \fig{varPicking}
rescale the numbers $[0-1]$ to $[0-15]$.
This range was chosen to include the values
that would lead to our desired behavior,
$[0-1]$, as well as many values that wouldn't.
The initial populations are chosen from 
$[.04-2.04]$.
Zero was disallowed because $(0,0)$ is a
trivial stable state in the model.

\fig{model} shows the model appropriately formulated
in the SPY language.
The change in the current populations of the two species
is calculated and returned to the main model executing code.
Notice that the two functions call the functions described
in \fig{varPicking}.

\fig{execution} shows the main control loop for the
competitive exclusion model.
It first sets the values that will be used for the
current simulation, lines 2-4.
Then it initializes the populations, line 6.
Lines 8, 9, and 12 write some useful information to disk.
Lines 15 and 16 call the two functions that comprise this
model, shown in \fig{model}.
Lines 17 and 18 adjust the current sizes of the two 
species.
Lines 19 and 20 ensure that model runs for at least 
{\tt REDO} time steps, but then stops when the populations
approach a steady state.
Finally line 21 is the worth function, which will be
described below.
This function is SPY's version of the objective function.


\subsubsection{Worth Functions}\label{sect:worth}
One practical detail of using SPY
that turned out to be more challenging than
expected was how to code the worth function.
The worth function is what SPY uses to decide
whether an instance belongs to the preferred class
or not, and hence plays the role of the objective function.
Tar4 needs discrete classes,
but since tar4 has a built-in discretizer,
worth functions can generate numeric class labels.
The worth function then must produce
a single numeric value that describes
how close to the preferred class any instance is.
The difficulty in formulating these functions
should be familiar to anyone who works
at an interface between natural languages and mathematics.
While it may be easy to describe a certain type of behavior
in a natural language, developing a mathematical formulation
can be difficult.

\paragraph{Competitive exclusion}
As mentioned above, in the dynamic population
model we were looking for a single behavior,
both species surviving at the equilibrium point.
So our worth function was
\begin{equation}
worth=(\text{if}~N_1 > 0 ~\text{then} ~1 ~\text{else} ~0) + 
(\text{if}~N_2 > 0 ~\text{then} ~1 ~\text{else} ~0)
\label{eqn:compexclusionworth}
\end{equation}
Note that the analytic analysis in \cite{biomodeling80}
shows if $N_1^0 > 0$ and $N_2^0 > 0$, one
species will always survive.
Equation \ref{eqn:compexclusionworth} therefore will only
return 1 or 2.

\paragraph{Animal neurons}
The worth function used in the animal neuron model
depended on which behavior we wanted SPY to look for.

\subparagraph{Phasic spiking}
A neuron exhibits phasic spiking when
its membrane potential ($v$ from equation \ref{eqn:neuron1})
spikes once and is then quiescent.
The model increments the {\it reset} counter every time
the reset condition is reached during each run.
The worth function used was
\begin{equation}
worth=\text{if} ~reset = 1 ~\text{then} ~1 ~\text{else} ~0
\end{equation}



\subparagraph{Tonic bursting}
A neuron exhibits tonic bursting when its membrane
potential spikes in repeated groups, where
each group has a very high frequency, but there
is a long quiescent period between each burst.
To classify this behavior the time between each
reset was recorded.
If the time since the last reset was below
a threshold, {\it low}, the counter {\it fast} was incremented.
If the time since the last reset was above
a higher threshold, {\it high}, (some times were therefore ignored
completely), the counter {\it slow} was incremented.
Finally at the end of the run the following worth
function was used
\begin{equation}
worth = \text{if} ~(fast > 10 ~\textbf{and} ~slow > 2) ~\text{then} ~1 ~\text{else} ~0
\end{equation}
The small number of {\it slow} resets delineated the groups, while
the large number of {\it fast} resets ensured that each group
had several resets.

\subparagraph{Phasic bursting}
Phasic bursting is similar to tonic bursting, except
there is only a single group of high frequency spikes.
Using the same counters as in the tonic bursting case,
the worth function was
\begin{align}
&worth = \notag \\
&\text{if} ~(fast > 4 ~\textbf{and} ~slow=0 ~\textbf{and} ~last\_spike > high) ~\text{then} ~1 ~\text{else} ~0
\end{align}
where {\it last\_spike} was the time since the last
reset when the simulation ended.
This worth function was particularly hard to develop.  
Recall that the definition of {\it slow} was based on 
the time between two resets.
With phasic bursting the quiescent period should last until
the end of the simulation, so there would be no 
reset to cause the {\it slow} counter to be incremented.
This last detail was difficult to formulate.

\subsection{Experimental Results}\label{sect:bioResults}

\subsubsection{Competitive Exclusion}
We ran 10 repeats with each repeat lasting 10 iterations.
Each iteration used only 100 instances.
The results of these experiments are summarized in
\ble{ceRelevant} and \ble{ceIrrelevant}.
\ble{ceRelevant} shows the bounds that SPY found on
the parameters that our analytic analysis showed
were relevant ($a_{12}$ and $a_{21}$).
\ble{ceIrrelevant} shows the bounds found on the irrelevant parameters
($\rho$, $u_1^0$, and $u_2^0$).
If there is a number in parentheses after the bound, 
it indicates the
number of times that particular bound was found.
The final column shows the average number of iterations
before SPY stopped adjusting the bounds on a particular parameter 
(lower or upper).
The original bounds on $a_{12}, a_{21}$ and $\rho$ were
$[0,15]$, while the original bounds on 
$u_1^0$ and $u_2^0$ were $[0.04,2.04]$
(0 could not be allowed because $(0,0)$ is a trivial stable
state).
\begin{table}
\center
\begin{tabular}{c|c|c|c|}
\hline
name & lower bounds & upper bounds & speed \\
$a_{12}$ & 0.0(10) & 0.2(6), 0.6(2), 0.8(2) & 1.8 \\
$a_{21}$ & 0.0(10) & 0.2(3), 0.3, 0.4(4), 0.8(2) & 1.2 \\
\hline
\end{tabular}
\caption{Relevant parameters in competitive exclusion model}\label{table:ceRelevant}
\end{table}

\ble{ceRelevant} shows that SPY was successful in finding ranges
that restrict the model to our desired behavior.  
Every repeat found restrictions that would always
lead to our preferred behavior (both species surviving 
in the stable state).
The speed to convergence was also very fast.
It took SPY only slightly more than 1 iteration to settle
on the range restrictions for $a_{21}$ and just under 2
iterations for $a_{12}$.
Analytically there is no difference between the influence
of $a_{12}$ and $a_{21}$, so we suspect the difference
in the speed of convergence for the two parameters
is an artifact of the way SPY utilizes the treatments
returned by tar4 or an artifact of the random search.
We can't call this experiment a complete success though, 
because the range restrictions were, for the most part,
too restrictive.
That is, the ranges SPY suggested eliminated a large
part of the parameter space that would lead to
the preferred class.

\begin{table}
\center
\begin{tabular}{c|c|c|c|}
\hline
name & lower bounds & upper bounds & speed \\
$\rho$ & 0.0(10) & 15.0(10) & 1.0 \\
$u_1^0$ & 0.04(10) & 2.04(10) & 1.0 \\
$u_2^0$ & 0.04(10) & 2.04(10) & 1.0 \\
\hline
\end{tabular}
\caption{Irrelevant parameters in competitive exclusion model}\label{table:ceIrrelevant}
\end{table}

\ble{ceIrrelevant} shows that SPY was completely successful
in ignoring irrelevant parameters.
Of the six irrelevant bounds SPY could have adjusted,
not a single one in ten trials was moved by SPY.

\subsubsection{Animal Neurons}
In this section we will discuss the results of the
experiments with the animal neuron model.
All experiments were run for 10 iterations, with
10 repeats for each behavior.
For this set of experiments a much larger batch size, 400,
was needed to allow SPY to converge to reasonable answers.

\begin{table}[t]
\center
\begin{tabular}{ccc}
parameter & lower bound & upper bound \\
$a$ & 0 & .1 \\
$b$ & 0 & 1 \\
$c$ & -80 & -40 \\
$d$ & 0 & 10
\end{tabular}
\caption{Original bounds on parameters for animal neuron model}
\label{table:neuronbounds}
\end{table}


There is a pair of figures for each behavior.
The first figure shows the average worth at the end of each iteration
for the ten iterations ran during the search
for all repeats that had a final average worth above .3
(in practice all repeats either had a final average worth 
close to 1 or close to 0).
(Recall a worth of 0 means that the 
behavior was not observed and a worth of 1 means that the 
behavior was observed.)
The second figure shows the range restrictions 
that SPY found by the last iteration for all
repeats that had a final average worth above .3.
Going from left to right, the four columns represent
the four independent variables in our model, $a,b,c$ and $d$.
Going from the bottom to the top, the different rows
represent the different successful repeats.
The top row shows the original bounds on the independent variables.
This makes it easy to see how much SPY restricted 
the range for each variable.
The original bounds of the variables are given in
\ble{neuronbounds}.
The same original bounds were used for all behaviors.
These bounds were set after reviewing
all the values used in \cite{izhikevich04} for the
20 different behaviors demonstrated.
The parameter settings used in 
\cite{izhikevich04} and \cite{izhikevich03}
for the different behaviors we are interested in
are listed in \ble{izhikevichsettings}.

\begin{table}[t]
\center
\begin{tabular}{lcccc}
& \multicolumn{4}{l}{parameter name} \\
behavior & $a$ & $b$ & $c$ & $d$ \\
\hline 
phasic spiking  & .02 & .25 & -65 & 6 \\
tonic bursting  & .02 & .20 & -50 & 2 \\
phasic bursting & .02 & .25 & -55 & .05 
\end{tabular}
\caption{Published parameter settings \cite{izhikevich04}}
\label{table:izhikevichsettings}
\end{table}


\paragraph{Phasic spiking}
Only four of the ten repeats found range
restrictions that exhibited phasic spiking.
The speed of convergence, when the search was
successful, was moderately fast.
\fig{PSworth} shows that for
the four successful searches, convergence happened
after 6-8 iterations.
A few qualitative statements can be made about
\fig{PSbounds}, which shows the final range
restrictions found by SPY.
The value of the $b$ parameter seems critical,
because all successful repeats found the same
narrow range.
The values of the $c$ and $d$ parameters seem
less critical, as a wide range (covering almost
the entire range in the case of parameter $d$)
of range restrictions was successful.
It is harder to say something about parameter $a$.
All successful repeats found very restrictive bounds
for $a$, but the bounds found do not coincide.
There may be a dependence between $a,c,$  and $d$
that SPY is finding, but this is only a conjecture.

\begin{figure}[hbt]
\center
\includegraphics[width=3.25in]{Figures/SPY_figs/neuronPS_worth.png}
\caption{Average worth for 10 repeated trials, trying to find phasic spiking}
\label{fig:PSworth}
\end{figure}

\begin{figure}[hbt]
\center
\includegraphics[width=6.0in]{Figures/SPY_figs/neuronPS2.png}
\caption{Final bounds for all successful repeats, trying to find phasic spiking}
\label{fig:PSbounds}
\end{figure}

\clearpage

\paragraph{Tonic bursting}
SPY was more successful at finding range restrictions
that exhibited tonic bursting.  Eight of the ten repeats
ended with an average worth very close to 1.
The speed of convergence was also faster while trying to
find tonic bursting.  It ranged from only 2 iterations to
5 iterations, as seen in \fig{TBworth}.
The larger number of successful repeats makes it easier
to make qualitative statements about the range restrictions
found, which are shown in \fig{TBbounds}.
First, we see something that was not observed with phasic
spiking: some of the bounds were not changed at all
by SPY.  For instance, many of the successful repeats
never lowered the upper bound on the $d$ parameter.
The same is true of the lower bound on the $c$
parameter and the upper bound on the $b$ parameter.
As with phasic spiking, the $d$ parameter has 
wide range restrictions, covering slightly more
than half the original range.
Tonic bursting is also less sensitive to the $b$ parameter 
than phasic spiking.
The noncoincidence of the found range restrictions on the 
$a$ and $b$ parameters again might be evidence of some
dependence between the parameters.

\begin{figure}
\center
\includegraphics[width=3.25in]{Figures/SPY_figs/neuronTB_worth.png}
\caption{Average worth for 10 repeated trials, trying to find tonic bursting}
\label{fig:TBworth}
\end{figure}

\begin{figure}
\center
\includegraphics[width=6.0in]{Figures/SPY_figs/neuronTB2.png}
\caption{Final bounds for all successful repeats, trying to find tonic bursting}
\label{fig:TBbounds}
\end{figure}

\clearpage

\paragraph{Phasic bursting}
SPY was more successful at finding good range restriction
for phasic bursting than it was at finding phasic spiking,
although not as successful as with tonic bursting.
Six of the ten repeats ended with an average worth close to 1.
The speed of convergence was slightly slower than
the speed for tonic bursting, 3-4 iterations, with one
outlier needing 8 iterations, as can be seen in \fig{PBworth}.
Phasic bursting appears to be very sensitive to
the $a$ parameter, with all successful repeats finding
the same narrow range, see \fig{PBbounds}.
The other three parameters show a wide range
of bounds, although it is difficult to say anything specific,
other than to note that the 2\nd successful trial has
the highest bounds on the parameters $b,c,$ and $d$.

\begin{figure}[hbt]
\center
\includegraphics[width=3.25in]{Figures/SPY_figs/neuronPB_worth.png}
\caption{Average worth for 10 repeated trials, trying to find phasic bursting}
\label{fig:PBworth}
\end{figure}

\begin{figure}[hbt]
\center
\includegraphics[width=6.0in]{Figures/SPY_figs/neuronPB2.png}
\caption{Final bounds for all successful repeats, trying to find phasic bursting}
\label{fig:PBbounds}
\end{figure}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "../thesis"
%%% End: 
