{'fraction_non_alphanumeric': 0.028101994724410807, 'fraction_numerical': 0.011637329103667052, 'mean_word_length': 4.418402615600187, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Applications of Discrete Mathematics for Understanding Dynamics of Synapses and Networks in Neuroscience

Mathematical modeling has broad applications in neuroscience whether modeling the dynamics of a single synapse or an entire network of neurons. In Part I, we model vesicle replenishment and release at the photoreceptor synapse to better understand how visual information is processed. In Part II, we explore a simple model of neural networks with the goal of discovering how network structure shapes the behavior of the network. To fully understand how visual information is processed requires an understanding of the way signals are transformed at the ribbon synapse of photoreceptor neurons. These synapses possess a ribbon-like structure capable of storing around 100 synaptic vesicles, allowing graded responses through the release of different numbers of vesicles in response to visual input. These responses depend critically on the ability of the ribbon to replenish itself as ribbon sites empty upon release. The rate of vesicle replenishment is thus an important factor in shaping neural coding in the retina. In collaboration with experimental neuroscientists we developed a mathematical model to describe the dynamics of vesicle release and replenishment at the ribbon synapse. To learn more about how network architecture shapes the dynamics of the network, we study a specific type of threshold-linear network that is constructed from a simple directed graph. The network construction guarantees that differences in dynamics arise solely from differences in the connectivity of the underlying graph. By design, the activity of these networks is bounded and there are no stable fixed points. Computational experiments show that most of these networks yield limit cycles where the neurons fire in sequence. We devised an algorithm to predict the sequence of firing using the structure of the underlying graph. Using the algorithm we classify all the networks of this type on five or fewer nodes.

To better understand vision, we wish to first understand how visual information is processed at the photoreceptor synapse. To do this we need to set up some necessary background information about the biology of these synapses. We will start by discussing neurons and how they communicate with each other. Then we will look more specifically at how photoreceptor neurons communicate. Finally we will describe the synaptic ribbon, a specialized structure in photoreceptor synapses, and discuss what is currently known about its role in the vesicle cycle and information processing.
Following the introduction, Chapters 2-4 discuss our work involving the ribbon synapse. These results were obtained in collaboration with experimental neuroscientists in the Thoreson Lab at the University of Nebraska Medical Center. Our contribution has been to develop theoretical models describing the dynamics of release and replenishment in the ribbon synapse. Results from Sections 2.1-2.3 were published in [36] and results from Sections 3.1-3.2 were published in [38].

Synaptic transmission
Neurons are cells involved in the transmission of information in the nervous system.
The neurons receive inputs from other neurons at the dendrites and once a threshold is reached the neuron can send a signal, often in the form of an action potential, down its axon to other cells. The pattern of action potentials codes the information being sent.
This information transfer between the two cells takes place at the synapse. When an action potential reaches the cell terminal, it triggers vesicle exocytosis.
Vesicle exocytosis is the process in which the vesicles dock and fuse with the cell membrane and release their neurotransmitters into the synaptic cleft, the space between the two cells. The area of the cell membrane where this occurs is referred to as the "active zone." The neurotransmitters then bind with the receptors on the postsynaptic cell, passing the information. For example, these receptors may open or close ion channels or activate second messenger systems.
Vesicles are recycled through endocytosis, which is the process by which vesicles are reformed using parts of the cell membrane and refilled with neurotransmitters.
These recycled vesicles then become part of the mobile vesicle pool inside the cell.

Photoreceptor neurons
We are particularly interested in studying synaptic transmission at photoreceptor neurons. Photoreceptor neurons are the first cells of the visual system. Photoreceptors are located in the retina and their function is to convert light into changes in membrane potential. Light is absorbed by membranous disks, located on the outer segments, containing photopigment. There are two main types of photoreceptors: rods and cones.
Rods are involved in night vision, motion detection, and peripheral vision and they are dense everywhere but the center of the eye. Cones are located in the center of the retina and are involved in color vision and detecting finer detail. Unless otherwise specified, all the experimental data in Part I refers only to cone photoreceptors, specifically in the aquatic tiger salamander.
In a conventional synapse the neuron responds to action potentials with discrete vesicle events. Photoreceptor cells instead respond directly to the absorption of photons by releasing vesicles constantly in darkness and slowing release as light increases, i.e. the cell is depolarized in darkness and an increase in light causes the cell to hyperpolarize. The graded responses given by photoreceptor cells allow for a quicker processing of information as well as a larger range of responses [19]. This graded release is facilitated by a structure called the synaptic ribbon, described in the next section.

Synaptic ribbon
The synaptic ribbons present in cone photoreceptors are plate-like rectangular 3 proteinaceous structures anchored to the inside of the cell membrane close to the Ca 2+ channels [33]. Daily and seasonal changes in the size, shape, number, and location of synaptic ribbons can occur based on light conditions [39]. Vesicle release at the active zone is controlled by the opening and closing of the calcium channels. In cones, less than three channel openings are required to cause the fusion of a single vesicle, which allows for precise timing of release to accurately reflect changes in light intensity [3].
The increase in intracellular Ca 2+ also speeds the replenishment of vesicles, allowing for sustained release [1].
Recall that in conventional synapses vesicles dock and fuse directly with the cell membrane. The vesicles in ribbon synapses are instead first collected on the synaptic ribbon. In the cone photoreceptors of the aquatic tiger salamander there are approximately 11 rows of 5 vesicles stacked on each side of the ribbon, for a total of 110 vesicles [2]. The vesicles become tethered to the ribbon via tiny filaments and then move along the ribbon towards the active zone. Not much is known about how the vesicles move down the ribbon to the active zone, but recent research posits that vesicles passively diffuse along the ribbon without an active transport mechanism 3 Synaptic ribbons in different cells may have different shapes. For example, ribbons in the auditory system can be spherical or ellipsoidal rather than rectangular ( [32], [17]).  [37]. [10]. Once the vesicles reach the bottom two rows of the ribbon they are considered part of the rapidly releasable pool (RRP). Experiments have shown that the ribbon may play a role in priming the vesicles for release [31]. As a result, the RRP can be released almost immediately following the opening of calcium channels. Once the RRP is depleted, additional vesicles from the reserve pool on the ribbon take their place. Empty sites on the ribbon are refilled by the mobile vesicles in the cell terminal. There are many theories regarding the function of the synaptic ribbon. The ribbon appears to support high rates of sustained vesicle release [33], but how the ribbon achieves this is still an open question. One theory is that the ribbon acts as a "conveyor belt" shuttling vesicles toward the active zone [23]. Another theory posits that it serves to hold the vesicles in contact with each other to facilitate multivesicular release via compound fusion [23]. Another theory asserts that the ribbon slows the delivery of vesicles, regulating the timing of release [16]. Yet another proposes that the ribbon functions to store the vesicles close to the active zone [45].

Questions about release and replenishment at the synaptic ribbon
With the goal of better understanding how visual information is processed at the photoreceptor synapse in mind, we ask some questions about the role of the synaptic ribbon in the vesicle cycle of this synapse.
The number of vesicles released is stimulus-dependent, with stronger stimuli resulting in more vesicles released. What causes this stimulus-dependence in the vesicle release? Does it depend solely on the probability of release or does it also depend the number of vesicles currently available on the ribbon? To answer this question we created a model of release and replenishment using experimentally measured quantities to predict the unknown quantities of pool size and release probability. This allowed us to independently predict the pool size and the release probability to determine which changes with the stimulus strength. See Chapter 2.
There is an upper limit on the rate of sustained vesicle release. What is the rate-limiting factor for release? Studies indicate that vesicle replenishment is the rate-limiting step in sustained release [16], so we take a closer look at replenishment.
Vesicles move randomly in the cell terminal without a directed movement toward the ribbon or active zone [24]. This may be due to the fact that ribbon synapses lack synapsins, proteins that help bind vesicles to the actin cytoskeleton, allowing the vesicles to diffuse freely [14]. Is this random motion of vesicles the rate-limiting step for replenishment? To answer this question we designed a three-dimensional random walk model of vesicle replenishment to calculate the replenishment timescale. We conclude that the random motion is not rate-limiting for replenishment. See Section 3.1.
It is known that Ca 2+ speeds the replenishment process [1]. By what mechanism does Ca 2+ speed replenishment? To answer this question we modified our random walk model to test whether Ca 2+ acted on the ribbon or on the vesicles. We conclude that Ca 2+ affects the probability of attachment at sites on the ribbon rather than directly affecting vesicles. See Section 3.2.
To further study replenishment we ask two additional replenishment-related questions: (1) How many vesicles collide with the ribbon per second? and (2) How long does it take for the ribbon to fill? We can use our random walk model of replenishment to answer both. See Section 3.3. Our random walk model does not take into account the geometry of the synaptic ribbon, so we designed a computational model to test the effects of ribbon geometry.
The computational model indicates that ribbon geometry does play a role in replen-ishment, so we also explore changes in local concentration near the ribbon and the effect of attachment probability on replenishment in an effort to explain the effects of geometry. The results appear in Section 4.

Chapter 2
Model of vesicle release and replenishment Photoreceptors respond to changes in light by releasing vesicles from the ribbon.
The amount of release depends on several key quantities: available pool size, release probability, and quantal amplitude. The available pool size, N , is the number of vesicles on the ribbon that are primed and ready for release. The release probability, P , is the probability that a vesicle on the ribbon will be released. The quantal amplitude, Q, is the postsynaptic influence of a single vesicle. We experimentally measure the response of photoreceptors to a given stimulus by measuring the postsynaptic currents 1 (PSCs) evoked in the postsynaptic cells onto which the photoreceptor synapses. Since the postsynaptic current is a linear sum of mini-EPSCs 2 [7], we can then estimate the number of vesicles released from the postsynaptic current using the quantal amplitude.
The amount of release depends on the stimulus, so which of N , P , and/or Q contribute to this stimulus-dependence? Changes in quantal amplitude occur on a longer timescale than our experiments [4,15,18,43]. Thus Q cannot change quickly enough to be the cause of the stimulus-dependent changes in postsynaptic response.
Stimulus-dependent changes in postsynaptic response are often due to Ca 2+ -dependent changes in P [36], but it is also possible that stronger stimuli allow Ca 2+ to spread further up the ribbon, effectively increasing N . With these possibilities in mind, are the stimulus-dependent changes then due only to changes in the release probability, P , or are they a result of changes in N as well?
In this section we will discuss a paradox that arises when asking this question. We then provide a model that estimates N and P independently, based on experimental data, allowing us to resolve this paradox. We will also describe a generalization of the model. The results in Sections 2.1-2.3 are published in [36]. The generalized model results in Section 2.4 are an unpublished extension of this work.

Dynamics of release and replenishment
Vesicles on the ribbon are released when the photoreceptor is stimulated and the amount released depends on the stimulus. As the vesicles are released, the empty sites on the ribbon are replenished by vesicles freely diffusing in the cell terminal. To study release and replenishment we model the available pool size during alternating periods of release and replenishment. Figure 2.1 shows a cartoon of the model. The variable A(t) tracks the number of vesicles on the ribbon at time t. When the stimulus is on, A(t) decreases due to vesicle release and when the stimulus is off, A(t) increases due to vesicle replenishment. When designing the model, we make several key assumptions.
We assume that no replenishment occurs during the periods of release. We also assume identical stimuli for the basic model, but the generalized model allows for multiple stimulus types in the same trial. By patching together the release and replenishment Release dynamics. The cumulative release at time t, c(t), is governed by the differential equation: where p s is the stimulus-specific probability of release, A i is the pool size at the beginning of the ith pulse, and τ r is the time constant of release. In salamander cones, the cumulative release curve can be fit by a two term exponential, one of the release time constants, τ r , is around 5 ms and the other is too long to be accurately measured in our experimental setup, so we omit it from the model. The timescale τ r regulates the release for strong stimuli (e.g. steps to -19 mV), and for weaker stimuli (e.g. steps to -39 mV) the time constant is made effectively slower by the release factor p s . Hence for strong stimuli we have p s = 1. Solving the differential equation yields Thus the release during the ith pulse is given by where P s = 1−e −ps∆t/τr . Note that for 25 ms steps to -19 mV, This is consistent with previous work showing that steps to -19 mV are strong enough to stimulate the release of nearly the entire pool of vesicles [2].

Replenishment dynamics
The cumulative replenishment at time t, a(t), is governed by the differential equation: where n is number of sites unoccupied at the end of a pulse and τ a is the time constant of replenishment. In salamander cones, replenishment is modeled with a two term exponential with time constants τ fast = 815 ms and τ slow = 13 s. Since the experiments occur on a much faster timescale than τ slow , we ignore τ slow in our model. Solving the differential equation yields Thus the amount of replenishment after the ith pulse is given by To study the release and replenishment dynamics we consider a pulse train experiment.
In this setup the cone is voltage-clamped near resting membrane potential and a steady train of pulses, i.e. voltage jumps, is applied to the presynaptic cell. Each pulse has duration ∆t ms, and the time between pulses is T ms. During a given pulse the voltage jumps up to a chosen voltage step (e.g. a step to -19 mV) and between steps the voltage returns to -79 mV (see Figure 2.2). The postsynaptic currents (PSCs) are measured in the postsynaptic horizontal cells. Since postsynaptic response is the result of a linear sum of independent quantal release events (mEPSCs) [7], we can use the postsynaptic measurements to estimate the number of vesicles released from the cones 3 . Once the release and replenishment reach an equilibrium where release is limited by replenishment, we can measure the limiting release. Our goal is to design a model that can predict maximum pool size and release probability using the first release and limiting release values measured during such an experiment.

An apparent paradox
Previous work uses a method of back-extrapolation to estimate the maximum pool size, A [25]. This method considers the cumulative release curve and fits a line to the steady state response that occurs when release is limited by replenishment. Backextrapolating to the time 0 gives an estimate of A. This method predicts that the maximum pool size is significantly smaller for weaker stimuli (see Figure 2.3). The amplitude of the releasable pool predicted by back-extrapolation is 80 pA for -39 mV, 105 pA for -29 mV, and 132 pA for -19 mV [36].
One of the pitfalls of the back-extrapolation method is that it assumes that the replenishment rate is constant. As we saw in Section 2.1, the replenishment rate is certainly not constant in salamander cones. This causes the method to underestimate  [36]. To predict maximum pool size we back-extrapolate from the steady state to t = 0.
the maximum available pool size since replenishment is faster when the ribbon has more available space. The method is close for the stronger stimuli because the pulse train stabilizes to the limiting release right away when exposed to a strong stimulus.
It also does not take into account the fast and slow replenishing sites. After the first pulse, slow sites don't have time to fill between pulses, so the limiting release reflects only the replenishment of the fast sites.

Using our model to predict pool size
In our model, we let both the maximum pool size A s and the release probability P s vary independently to determine which causes the voltage-dependent changes in vesicle release. We want to estimate both A s and P s in terms of the measured first release (R 1 ) s and limiting release R s for each stimulus s.

Derivation of pool size and release probability formulas
Let A i be the pool size at the beginning of the ith pulse, c i (t) be the cumulative release t milliseconds into the ith pulse, and a i (t) be the cumulative replenishment t seconds into the ith pulse. Then we can compute A i by taking the pool size at the beginning of the previous pulse, subtracting the release during that pulse, and adding the amount replenished before the ith pulse, i.e.  [36]: A(t) keeps track of pool size at time t. Stimulus pulses of duration ∆t cause vesicle release and between pulses we have replenishment periods of duration T . Note that the pool size decreases during the pulses due to vesicle release and increases in between the pulses due to replenishment. The maximum possible pool size is denoted by A.
Solving the recursion we get The details of solving the recursion appear in Lemma 2.1. Taking the limit as i → ∞ gives which represents the pool size at the beginning of each pulse during the steady state.
Since R i = P s A i is the amount released during the ith pulse, then the limiting release is given by Solving for A s yields Also, note that (R 1 ) s = A 1 P s = A s P s . Solving for P s and substituting into our equation for A s gives . (2.12) Using the fact that (R 1 ) s = A s P s , we can also find an expression for P s ,

Generalized pulse trains
We consider a generalized pulse train with n periods of release and replenishment. The length of the th period is ∆t . We do not necessarily alternate between release and  replenishment. Once we have cycled through all n periods we started at the beginning and repeat the periods in the same order. We denote the available pool size at the beginning of the th period of the ith cycle by A i, and assume the ribbon is full at t = 0, so A = A 1,1 is the maximum available pool size. See Figure 2.5. Let a i, (t) be the total change in pool size t seconds into period of cycle i with a i, (0) = 0. If index j denotes a release period then the dynamics is governed by If index k denotes a replenishment period then the dynamics is governed by

Setting up and solving the recursion
Note that this setup gives Then we get In release periods a i,j (t) = −A i,j (1 − e −t/τ j ) and in replenishment periods a i,k (t) = Let α = e −∆t /τ so we can rewrite the terms in the sum using Then ordinary generating functions gives a solution of Note that since b = n =1 α < 1, then Thus we have a limit cycle.
For s.t. θ( ) = 0, we have that the total release during that period is given by Lemma 2.2 gives a closed formula for A i, in terms of A i,1 , so

Special cases of the generalized model
Release only. Consider the case where each period is a release period. Then L = ∅.

So for the recursion we have
which has solution When we take the limit we get Replenishment only. Consider the case where each period is a replenishment period. Then L = {1, . . . , n}. So for the recursion we have which has solution When we take the limit we get Alternating release and replenishment periods. Consider the case where we have alternating periods of release and replenishment, starting with release. Then L = {2, 4, . . . , n}. So for the recursion we have which has solution When we take the limit we get We can calculate the total release during release period by Alternating release and replenishment, original model. In the case of our experiment we have α 1 = α, α 2 = β, L = {2}, and n = 2. So for the recursion we have which has solution Taking the limit we get , which matches our prediction from Section 2.3.1 (without the f ).
Two different release periods, but identical replenishment. Here we have alternating periods of release and replenishment with two different release periods, starting with release, so L = {2, 4}. We also assume that all of the replenishment periods have the same dynamics. Let β := α 2 = α 4 . So for the recursion we have which has solution Taking the limit we get

Supporting lemmas
In this section, we give the proof of the recursion in Equation 2.17. We also give proofs of the supporting lemmas used to prove Proposition 2.4 (Equation 2.19).

Lemma 2.1. The solution to the recursion
In the following proof we use a standard generating function technique for solving recursions found in [42].
Proof. To solve using generating functions we first multiply the recursion by x i and sum over i ≥ 2 to get: Solving for A(x) yields: Proof. We induct on . Note that for = 1 we have So the result holds for = 1. Let > 1 and assume that the result holds for all smaller . Then by the induction hypothesis Proof. We know that Then we have Chapter 3

Random walk model of vesicle replenishment
Previous work indicates that vesicle replenishment is the rate-limiting step in sustained vesicle release [16], so in this chapter we take a closer look at the replenishment process.
What part of the replenishment process limits release? In Section 3.1, we design a simple random walk model to theoretically predict the time constant of replenishment, τ a , initially discussed in Section 2.1. Using the model we can determine which parameters affect replenishment. We discover that τ a relies on four fundamental parameters: vesicle diffusion, vesicle concentration, vesicle size, and the probability of attachment to the ribbon. The model predicts an exponential replenishment curve with time constant where D is the vesicle diffusion coefficient, ρ is the vesicle concentration, δ is the diameter of a single vesicle, and s is the attachment probability. The nature of vesicle movement within the synapse leads us to ask if the random diffusion of vesicles is rate-limiting. We compare experimental data with our model results and conclude that diffusion is not, in fact, rate-limiting.
Further exploring replenishment in Section 3.2, we introduce two variations of the original model to investigate the role of Ca 2+ on replenishment. The results in Sections 3.1 and 3.2 are published in [38].
We can also use the model to calculate several other quantities of interest: how many vesicles hit the ribbon per second (hit rate) and how long it takes to fill up the ribbon (expected waiting time). The derivations in Section 3.3 are unpublished.

Replenishment timescale
In this section we wish to answer the question: Is vesicle diffusion a rate-limiting step for replenishment? We first discuss how to measure the replenishment curve experimentally. Then since vesicles move randomly in the cell terminal without a directed movement toward the ribbon or active zone [24], we create a random walk model of vesicle movement and replenishment.

Paired pulse recordings
To experimentally measure the replenishment of vesicles onto the synaptic ribbon, we use paired pulse recordings. First, when the ribbon is full, a large voltage jump, or pulse, is applied to the cell (similar to the pulse trains in Section 2.2.1) and vesicle release is measured. The size of the voltage jump is such that all vesicles from the ribbon are released. After t seconds a second pulse of the same amplitude is applied and vesicle release is again measured (see Figure 3.1). This is repeated for multiple values of t to approximate the replenishment curve. For each inter-pulse interval, t, we plot the ratio R 2 /R 1 where R 1 is the release on the first pulse and R 2 is the release on the second pulse (see Figure 3.1). The ratio, R 2 /R 1 , can be thought of as the percentage of vesicles replenished. In the case of salamader cones, the replenishment curve can be fit with a double exponential with time constants τ fast =815 ms (76%) and τ slow =13 s [38]. Figure 3.1: Panel A shows an example of two paired pulse recordings with interpulse intervals of 500ms and 2s. Note that the longer interpulse interval gives more time for the ribbon to replenish and hence the the second pulse is larger in the 2s trial. Panel B shows the replenishment curve. The horizontal axis gives the interpulse interval and the vertical axis gives the ratio of the two responses. This ratio can be thought of as the percentage replenished. Adapted from [38].

A B
In the next section our goal is to predict the time constant of replenishment theoretically using a random walk model.

Derivation of the replenishment time constant
To answer this question, we developed a three-dimensional random walk model. We modeled the vesicle motion in the synapse by spherical vesicles undergoing random walks on a rectangular lattice of spacing δ. During each time step, ∆t, every vesicle moves to an adjacent lattice site in each dimension. We update each of the three dimensions simultaneously, resulting in a diagonal move overall. The (macroscopic) diffusion coefficient, relates δ and ∆t in the (microscopic) random walk model, so these quantities cannot be chosen independently [6]. Moreover, we would like to assume that each lattice site can be occupied by at most one vesicle, and that the occupation probabilities for distinct lattice sites are independent. These two assumptions can only be satisfied if we choose δ to be equal to the vesicle diameter. We use p to denote the probability that a given lattice site (or tethering site) on the ribbon will become occupied in a given time step. If we assume the vesicles are distributed randomly and uniformly within the cell, the probability of a given lattice site being occupied is independent from one time step to the next. For a lattice site far from the ribbon, this probability is simply given by the vesicle density per lattice site, ρδ 3 , where ρ is the overall density of vesicles inside the cell. Since the ribbon sites can only be accessed from one side, for ribbon sites we have collision probability Let s be the attachment probability, the probability that a vesicle that comes into contact with the ribbon will "stick." 1 Then sp is the probability of a vesicle actually sticking to a ribbon site in a given time step. Thus the probability of having to wait at least t seconds before a ribbon site is "permanently" occupied is: with t/∆t giving the total number of time steps that have elapsed in t seconds. Note that 1 − sp is the probability that a given lattice site on the ribbon is not occupied permanently in a given time step.
Now we make a crucial approximation for P (t), which is valid for sp << 1. 2 The approximation stems from the fact that ln(1 + x) ≈ x for |x|<< 1. To use it, we first take the natural log of the P (t) equation, and then plug in ∆t = δ 2 /2D and p = ρδ 3 /2: Exponentiating both sides we obtain 1 i.e. become tethered to the ribbon until release, not drifting away at a future time step. 2 In fact, the approximation is still quite good up to values of sp ∼ 0.1.
Solving for P (t) without making the approximation we get Next, observe that the expected number of ribbon sites that are filled at time t, assuming all sites are empty at t = 0, is given by The second equality is obtained using a familiar variant of the Binomial Theorem.
Recall that n m=0 n m x m y n−m = (x + y) n , by the Binomial Theorem. Differentiating with respect to x: Now, letting x = 1 − P (t) and y = P (t), we obtain Finally, multiplying both sides by 1 − P (t) we obtain the desired result. Note that each term in the sum corresponds to the probability that exactly m sites are "permanently" occupied at time t, weighted by m. On the other hand, given that each of the n ribbon sites has an occupation probability of 1 − P (t) at time t, it is intuitive that the expected number of occupied sites at this time is a(t) = n(1 − P (t)).
Using the approximate expression for P (t) we derived above, we obtain, a(t), the expected number of vesicles on the ribbon at time t, in terms of our fundamental constants:

Comparison of model predictions with data
We wish to determine whether diffusion is rate-limiting for replenishment and thus sustained release. We answer this question by comparing the experimentally measured time constant with the model predictions. Recall that in Section 3.1.1, we fit an exponential replenishment curve with two time constants to the data with τ fast =815 ms (76%) and τ slow =13 s.

Role of calcium in replenishment
In this section we use variations on our model to test two different mechanisms by which calcium (Ca 2+ ) and calmodulin (CaM), a calcium-binding messenger protein, might govern the attachment probability, s. It is known that Ca 2+ speeds replenishment [1]. However, the mechanism by which this occurs is unknown. Data suggests that Ca 2+ /CaM do not accelerate vesicles from the top of the ribbon to the release sites, nor do they increase the fusion rate at the membrane [38]. Increased intracellular Ca 2+ does not affect the mobility of vesicles in the terminal [24,14]. Hence D and ρ would not be affected by calcium. It appears that vesicle size (quantal amplitude) is also not affected by calcium [38], so we posit that Ca 2+ /CaM increases the attachment probability, s. In this section, we will use two variations of the model to test two hypotheses regarding the role of Ca 2+ /CaM. The first variant, which we call Model 1, assumes that Ca 2+ /CaM acts as a switch making some vesicles more "sticky" than others. In the second variant, Model 2, we assume Ca 2+ /CaM again acts as a switch, but this time on the ribbon tethering sites, making some ribbon sites more "sticky" than others but leaving the vesicles unchanged. Perhaps surprisingly, these two models

Calcium affects vesicles
The rest of the model remains unchanged. In particular, we still have  Since attachment probabilities are different for different ribbon sites, we must use different expressions for P (t):

Calcium affects the ribbon
The result is that expected number of vesicles on the ribbon at time t is given by the sum of two terms:

Comparison to experimental results
Recall that paired pulse experiments showed that the replenishment curve can be fit with a double exponential. Table 3.3 shows the results of experimentally decreasing Ca 2+ /CaM using BAPTA, nifedipine, Calmidazolium, and MLCK [38]. When the fast timescale is constrained to match the control, we see that the fit is comparable to the unconstrained case, but the percentage of fast-replenishing sites is much lower.
Thus, inhibiting Ca 2+ /CaM in these experiments caused slight changes in the fast timescale, but more substantial changes to the amplitude of the fast component. This is consistent with the predictions seen in Model 2, where we have two time constants and inhibition of Ca 2+ /CaM causes the amplitude of the fast component to decrease.
We conclude that Ca 2+ /CaM more likely acts on ribbon sites rather than vesicles. It is possible that vesicles are affected by Ca 2+ /CaM as well, but changes in the fast timescale were not consistent across trials.

Hit rate
In this section we compute the hit rate, i.e. the number of vesicles coming in contact with the ribbon per second. One way to compute the hit rate is to do so macroscopically by first calculating the flux near the ribbon and then multiplying by the surface area of the ribbon. We consider the concentration of vesicles to be zero on the ribbon.
The concentration of mobile vesicles not attached to the ribbon is ρ as before. Thus over the distance δ (one lattice step from a ribbon site to a nonribbon site) we have a change in concentration from ρ to 0 giving us a flux of

Expected waiting time
To calculate the expected waiting time, T wait , to fill all n lattice sites on the ribbon, we first consider an individual lattice site. Let P (t) be the probability that we wait at least t seconds to fill the given lattice site. Then 1 − P (t) is the probability that the given lattice site fills before t seconds have passed and r(t) = (1 − P (t)) n is the probability that all n sites have filled before t seconds have passed. Hence the probability we wait exactly t seconds is r (t)dt.
Thus the expected waiting time T wait is given by where τ a is the vesicle accretion timescale and H n = n k=1 1 k is the nth harmonic number.
This result is proven in Lemma 3.1. Proof. We have r (t)dt = n τa e −t/τa (1 − e −t/τa ) n−1 dt, so Using the identity (k + 1) n k + 1 = n n − 1 k and reindexing, we can rewrite τ a n−1 Then, letting u = 1 − x, we have, Thus,  Expected waiting time in Model 2 Recall that in Model 2, the stickiness occurs in the ribbon sites instead. We have n A ribbon sites with attachment probability s A and n B sites with attachment probability s B . We know that the probability of having to wait at least t seconds before a ribbon site in population A is occupied is given by Dρδs A and the probability of having to wait at least t seconds before a ribbon site in population B is occupied is given by is the probability that all n = n A + n B sites have filled before t seconds have passed.
Hence the probability we wait exactly t seconds is r (t)dt. Now, we have that the expected waiting time for the ribbon to fill is

Chapter 4 Computational model
The

Effect of ribbon geometry on replenishment
Since the theoretical model does not take into account the geometry of the ribbon we use our computational model to approximate the replenishment curve in the case where we have a rectangular ribbon attached to the edge of the cell space. Recall that in the analytical model we have a factor of 1/2 that represents the fact that the ribbon sites are only accessible from one side. When using the analytical model to predict replenishment in nonribbon cases, we leave out the factor of 1/2 since these sites are accessible from all sides. In the case of high attachment probability and low concentration, we hypothesize that once the ribbon starts to fill up, the local concentration near the ribbon decreases causing the ribbon to fill slower than predicted. To test this we calculate the concentration of vesicles close to the ribbon and far away from the ribbon.  Movement on ribbon and vesicle fusion. The random walk model does not take into account the movement of vesicles along the ribbon. As more becomes known about this process, we would like to incorporate this step into the model. This model also does not take into account vesicle release. Adding these features will allow us to explore more questions regarding the function of the ribbon.   Figure 4.2) are averaged over 100 trials and the gray area represents one standard deviation from the mean. We show trials for a low (300 vesicles/ µm 3 ) and a high concentration (2300 vesicles/ µm 3 ) as well as three different attachment probabilities (0.1, 0.5, and 1). Note that the theoretical and computational models appear to closely match across a variety of parameters, as expected in the nonribbon case.   Figure 4.2) are averaged over 100 trials and the gray area represents one standard deviation from the mean. We show trials for a low (300 vesicles/ µm 3 ) and a high concentration (2300 vesicles/ µm 3 ) as well as three different attachment probabilities (0.1, 0.5, and 1). Note that the computational model shows the greatest deviation from the theoretical prediction in the low concentration/high attachment probability case.

Part II Neural Sequences in
Threshold-Linear Networks

Introduction to Part II
Part II focuses on neural networks and the interplay between network connectivity and neural activity. In particular, we are interested in studying how network structure shapes the behavior of the network.
To do this, we study the dynamics of a combinatorial family of competitive threshold-linear networks constructed from simple directed graphs (the CTLN model) as defined in [21]. This family of networks is particularly well suited for our study because the network construction guarantees that differences in dynamics arise solely from differences in the connectivity of the underlying graph. This allows us to focus on the properties of the graphs themselves when trying to predict the behavior of the corresponding network. This robust family of dynamical systems exhibits several different nonlinear behaviors including limit cycles, quasiperiodic attractors, and chaos. In this part, we begin by giving some background about competitive thresholdlinear graphs and the CTLN model. We then use the CTLN model to study how the graph structure affects the resulting dynamics. Computational experiments show that most CTLN networks yield limit cycles. We present an algorithm that uses the structure of the underlying graph to predict the sequence of firing of neurons in the limit cycle. Our algorithm predicts the sequence correctly for most small graphs, but sometimes fails for certain classes of larger graphs. To gain further insight into how the structure of the underlying graph shapes the dynamics, we classify the behavior we see for small networks (n ≤ 5 nodes) arising from oriented graphs. Both of these results work towards the larger goal of better understanding high-dimensional nonlinear dynamics.

Threshold-linear networks
The CTLN model is a specific type of threshold-linear network. Neuroscientists use threshold-linear networks to model recurrent neural networks [27]. These networks are thought to be involved in perception and memory processes [11]. Development of the mathematical theory behind threshold-linear networks is ongoing [44,11,21].
Definition 5.1. A threshold-linear network on n neurons is defined by the following system of differential equations: where x i is the firing rate of the ith neuron, W is the matrix of connection strengths, θ ∈ R is the external drive to the network, and [y] + = max{0, y} is the threshold nonlinearity.
In our neural network context, the −x i represents the leak term and guarantees that the activity of neuron i will die out in the absence of other inputs. Inside the nonlinearity we have a sum of inputs from all other neurons weighted by the connection strengths. In inhibitory networks (W ij ≤ 0), the parameter θ must be positive in order for the nonlinear term to be nonzero.
We study the behavior of threshold-linear networks of n neurons as they are one of the simplest examples of a nonlinear system of ordinary differential equations. In particular, we are interested in studying the dynamics of a competitive threshold-linear network defined from a simple directed graph. with the added restriction that W ij ≤ 0 and W ii = 0 for all i, j = 1, . . . , n and θ > 0.
In the next section we will describe the CTLN model, which is a particular type of competitive threshold-linear network.

Description of the CTLN model
The Combinatorial Threshold-Linear Network model (CTLN model) was first introduced by Curto et al. in [8] and further explored by Morrison et al. in [21]. This model was designed as a way to study high-dimensional nonlinear dynamics without using a linear approximation [21]. Linear models are limited as tools for approximation as they do not demonstrate complex behaviors such as limit cycles, multistability, and chaos. The nonlinearity in the CTLN model captures the full range of nonlinear behaviors, but is still simple enough that it is possible to develop a corresponding mathematical theory. In this chapter we describe the CTLN model and necessary background. with only two values for the inhibitory connection strengths. For any δ > 0 and 0 < ε < 1, the n × n connectivity matrix W is given by where i ← j represents a directed edge from node j to node i in the graph G and i j means that no such edge exists in G [21].
The CTLN networks are therefore fully inhibitory, with each node acting on its neighbors by quieting their activity. A biological motivation for the model is shown in Figure 5.2. Inhibitory interneurons (gray circles) inhibit all neighboring excitatory pyramidal cells (colored triangles) equally [8]. Connections between excitatory neurons therefore have two strengths. A directed edge represents an overall connection strength of −1 − ε, i.e. inhibition has been weakened by an excitatory connection. Lack of an edge represents an overall connection strength of −1 − δ. The following theorem from [21] gives constraints on the graph and the relationship between δ and ε that guarantees bounded activity, but disallows stable fixed points.
A B If ε < δ 1 + δ , then the network has bounded activity and no stable fixed points.
By forbidding stable fixed points, Theorem 5.4 guarantees that the activity of the network is either oscillatory or chaotic. Computational experiments show that most of the small networks for which Theorem 5.4 holds exhibit limit cycles where the neurons often appear to fire in sequence. The construction of the CTLN model guarantees that any differences in dynamics arise solely from differences in the underlying graph.
Our goal is to use the structure of the graphs to predict the resulting sequences. In the following sections we will assume θ = 1, ε = 0.25, δ = 0.5, and we will focus on oriented graphs with no sinks. Theorem 5.4 holds in these cases. The next section shows some examples of these networks and their behaviors.

Examples and behaviors
In this section, we will explore examples of CTLN networks with a small number of nodes. We will start with the simplest example of a network satisfying the conditions of Theorem 5.4, a three-cycle as seen in Figure 5   Note that it is possible for a network to have multiple limit cycles or chaotic attractors.
For example, the network in Panel D has four different chaotic attractors, only one of which is shown.
Could we have predicted these dynamics by looking at the graphs? The limit cycles in panels A and B, each follow a cycle in the corresponding graph. This is common in smaller graphs, but for larger n we have seen examples where this is not the case. Also, what happens when there is more than one 5-cycle in the graph? Note that the graph in Panel B has two 5-cycles, 12534 and 15423, but the network has only one limit cycle. Why does the network preferentially choose one 5-cycle over the other? Figure 5.5 shows an example on n = 7 nodes. Note that the limit cycle shown in the activity trace corresponds to a 6-cycle in the graph. Why does node 2 stop firing? It receives input from three other nodes while nodes 3 and 5 only receive two inputs each. Additionally, there are multiple 7-cycles in the graph. Why doesn't the limit cycle correspond to one of these 7-cycles? These are all questions that our algorithm must address. The algorithm must be able to discern between limit cycles and chaos, predict the firing sequence for limit cycles, and must also identify which nodes stop firing, if any. show some examples of oriented graphs on n = 5 nodes. Panels A2-D2 show the dynamics of the networks corresponding to the graphs. Panel A2 shows a typical limit cycle. Note that nodes 4 and 5 fire at a much lower rate than nodes 1, 2, and 3. Panel B1 shows an example of a balanced subgraph. Note that each node has indegree 2 and outdegree 2. Panels C1-C2 show an example with synchronous firing: nodes 1,4, and 5 fire at the same rate. This is caused by the graph automorphism in C1. Panels D1-D2 show an example of a network with a chaotic attractor. Note that we show a longer trial in D2 to show the chaotic behavior. The limit cycle has a sequence of 634517 as indicated above the dynamics. Note that node 2 (red) stops firing. Adapted from [21].

Chapter 6
Sequence prediction algorithm In this section, we will discuss an algorithm we designed to predict the neural sequence from the graph for CTLN networks. The basic premise of this algorithm comes from the idea of removing the "weakest" node and looking at the dynamics of the remaining network. Once we know how the smaller network behaves we work to figure out a way to tell where the deleted node fits in the sequence. We start with a description of the algorithm, discuss some examples, and then state some conjectures about when the algorithm is successful. We will look first at the case of tournaments 1 without sinks and then examine oriented graphs without sinks.

Description of the algorithm
The algorithm has two separate phases. In the deconstruction phase, we will first deconstruct the graph by deleting one vertex at a time, keeping track of the order of deletion. Then in the reconstruction phase, we start with a base sequence based on our deconstruction and rebuild the neural sequence by adding back in the deleted Deconstruction phase. At each step of the algorithm we delete one of the vertices of G with smallest indegree such that the resulting reduced tournament has no sinks.
We continue to delete vertices until we can no longer do so. This occurs when the resulting tournament is a three-cycle (see Proposition 6.1). We refer to this three-cycle as the core cycle. At each step we record the current tournament and the vertex we deleted. See Figure 6.1 for an example. Note that the choice of vertex to delete is not necessarily unique, so it is possible for the algorithm to output multiple sequences.
Reconstruction phase. To reconstruct the sequence, we start with the three-cycle from the last step of the deconstruction phase. Recording the vertices in order of the three-cycle we insert the other vertices into this sequence in reverse order of deletion.
Proceeding backwards through the list of deleted vertices, we add each vertex back into the sequence following the vertex that feeds into it in the graph from the preceding step in the deconstruction. If more than one vertex feeds into the vertex to be added, we look at the subgraph induced by these possibilities and if one of these possibilities is a sink in the induced subgraph, we place the vertex to be added after the sink in the sequence. See Figure 6.1 for an example. Note that it is possible for the algorithm to fail if there are two or more edges feeding into the node we are adding back in.
Node death. When reconstructing the sequence, there are rules to predict the death of a node. If the vertex we delete at a given step has indegree zero then we do not add that node back in during reconstruction. If the vertex we delete at a given step has indegree one with the one edge coming from a vertex not in the core cycle then we do not add that node back in during reconstruction. of edges can also be given by m 2 since it is a tournament, so m = m 2 . Solving for m gives m = 3. Since we have outdegree 1 at every node, the current graph is a 3-cycle. Also note that if the graph at the current step is a three-cycle, then deleting any vertex, will result in a graph with two vertices and a directed edge between them.
Thus one of the vertices is a sink, so the deconstruction phase terminates at at step n − 3.
6.2 Performance of the algorithm Proposition 6.2. For ε = 0.5, δ = 0.25, and θ = 1, the algorithm correctly predicts the neural sequence for tournaments without sinks on n ≤ 5 nodes.
The proof of this proposition is done by checking each case computationally. We also note that for oriented graphs, the algorithm appears to predict which neurons will have high firing rates. For small examples, we often see three neurons with higher firing rates than the remaining neurons. These three high-firing neurons appear to correspond to the neurons in the core cycle predicted by the algorithm. Further exploration of tournaments lacking sinks on n > 5 nodes indicates that there are networks where the algorithm does not correctly predict the behavior, often in the form of spurious predictions or incorrectly predicting neuron death. Analysis of these graphs indicates that the networks for which the algorithm fails appear to have the common property of having a balanced subgraph on n ≥ 5 nodes or are an outerneuron construction.  We then add a directed edge from the pseudo-sink to the pseudo-source to guarantee that the new graph on n + 2 vertices has no sinks.
We have looked at all tournaments having no sinks on up to n = 7 vertices. Proposition 6.2 gives that the algorithm works for the 11 such graphs on n ≤ 5 vertices. On n = 6 nodes there are 44 graphs and on n = 7 nodes there are 400 such graphs. The algorithm fails for only two of the n = 6 graphs, one with a spurious prediction and one with a spurious deletion. Both of these graphs have a balanced subgraph on 5 vertices. For n = 7, if we ignore graphs with an outerneuron construction and graphs with balanced subgraphs with n ≥ 5 vertices (153 total), the algorithm only fails for 5 out of the remaining 247 graphs. Additionally, of the 153 graphs having an outerneuron construction or a balanced subgraph with n ≥ 5, 62 of these are still correctly predicted by the algorithm.

Extending the algorithm to oriented graphs
We have also explored the success of the algorithm on oriented graphs without sinks on n ≤ 5 nodes. Recall that unlike tournaments, oriented graphs do not require an edge between every pair of vertices. As a result we are not guaranteed that the algorithm will terminate at a three-cycle for oriented graphs. In fact, we have seen examples on n = 5 where the algorithm terminates in a 4-or 5-cycle. A comprehensive list of oriented graphs without sinks on n ≤ 5 nodes appears in Appendix C. The algorithm correctly predicts the behavior in all but 6 of the 160 total networks on n ≤ 5 nodes. Using the algorithm on the Graph #147, 148, 149, 152, and 158 predicts the correct sequence, but also produces a spurious prediction. For example, the network corresponding to Graph #147 has a limit cycle with sequence 12 (45)3 where the parentheses indicate that neurons 4 and 5 fire synchronously. The algorithm makes three predictions: 12453, 12543, and 12534. The first two predictions result in a correct final sequence of 12(45)3, but 12534 is a spurious prediction. Using the algorithm on Graph #153, a node is deleted that does not die, which keeps the algorithm from predicting the synchronous activity.
One weakness that arises when using the algorithm on oriented graphs rather than tournaments is that it is possible to disconnect the graph during deconstruction. Since the reconstruction rules will not necessarily make sense in this case, additional rules will be necessary for the deconstruction of oriented graphs to avoid breaking the graph into multiple components. We also still sometimes make incorrect predictions in the case of graphs with the outerneuron construction and/or a balanced subgraph. The next section shows a comprehensive study of all the oriented graphs on n ≤ 5 vertices in order to investigate ways to adjust the algorithm for oriented graphs.
6.4 An application: classification of oriented graphs on n ≤ 5 To investigate why the algorithm fails in some cases, we perform an exhaustive study and classification of oriented graphs on n ≤ 5 nodes. See Appendix C for a complete list of graphs and their classification. If we look at the graphs for which the algorithm failed, we see that they all fall into a category where there are two or more different n = 4 subgraphs possible in the first step in the algorithm. Using this classification we hope to start refining the algorithm to work in more generality.
In the next section we also sort the graphs by dynamics. Note that the graphs in the same entry in the dictionary often also appear in the same category of the classification in Appendix C.

Dictionary of attractors for n ≤ 5
In this chapter we create a dictionary of graph behaviors. We sort the graphs into groups based on the dynamics of the corresponding network. Each entry in the dictionary corresponds to a specific limit cycle or chaotic attractor. AT denotes "attractor type." For each entry, we show a representative graph for that particular attractor type along with its dynamics. To make the dynamics plots we use ε = 0.5, δ = 0.25, and θ = 1 as before. Initial conditions used are listed above the dynamics plots. We use an asterisk (*) to indicate which initial condition we used to create the plot. The sequence listed is for the representative graph only. We underline low-firing neurons and synchronous neurons are in parentheses. All graphs listed in that entry exhibit the same behavior (up to permutation) as the representative graph for some initial condition. If a graph has more than one attractor we annotate the graph number with ic1, ic2, etc. The labelling of the graphs corresponds to the numbering in the catalogue found in Appendix C. The compilation of the dictionary was carried out in collaboration with Katherine Morrison.

AT-18 (chaotic)
Rep. graph 150. All graphs: 150_ic1, 150_ic2*, 150_ic3, 150_ic4. Balanced subgraphs and outerneuron constructions. For the subset of graphs whose sequences cannot be predicted by our current algorithm, we need to adjust the current algorithm or design a new algorithm to handle these cases.
Refining the algorithm for oriented graphs. Using the classification of oriented graphs on n ≤ 5 nodes, we can look for commonalities in the graphs where the algorithm fails and adjust the algorithm accordingly.
Exploring larger networks. After classifying the networks and behaviors for oriented graphs on n ≤ 5 one of the next steps is to take a closer look at graphs with n > 5.
Proving conjectures about algorithm. We have many conjectures about the prediction algorithm, so we would like to look for proofs or counterexamples. Conjecture 6.5. In an oriented graph with no sinks on n = 5 nodes, two nodes fire synchronously if the algorithm predicts two sequences that are identical except for with the two nodes switched and neither node is part of the core cycle. Conjecture 6.6. During the reconstruction phase, if node u has indegree zero in the corresponding subgraph, then this node dies. Conjecture 6.7. During the reconstruction phase, if node u, with indegree 1 in the corresponding subgraph, is added to the sequence following a node that was not part of the base sequence then u dies, unless u is part of a balanced subgraph.

Appendix A
Code for random walk model of replenishment This appendix gives the Matlab code for the random walk model of replenishment.

• ReconstructCycle
The main function, run Algorithm takes the transposed adjacency matrix of a graph as an input and outputs the list of predicted limit cycles for the corresponding network. This appendix includes a list of all the oriented graphs on n ≤ 5 vertices with no sinks.
The graphs on n = 5 vertices are sorted by which n = 4 subgraph(s) they reduce to in the first step of the algorithm from Chapter 6. We group the n = 4 graphs into four classes: #2, 3, and 4 reduce to a three-cycle, #5 is a 4-cycle, #6 and 7 have three strong neurons and one weak, and #8 has two distinct limit cycles. For graphs whose networks exhibit more than one behavior, like #8, we list the initial conditions for each attractor type.  The left column gives the indices of the possible n = 4 subgraph(s) that appear after step one of the algorithm, the middle column shows the indices of the graphs that have this reduction (based on the classification indexing), and the right column shows which attractors these graphs have (based on the dictionary). Subgraphs 2, 3, and 4 are grouped together since they have the same behavior and subgraphs 6 and 7 are similarly grouped.