Theory One hypothesis on how ideas are structured in the brain is that groups of cells form a unit called a cell assembly. The cells within the assembly are connected with strong excitatory synapses. Thus, when enough activity occurs in the neurons within the assembly, then the entire assembly will have the ability to activate. An active assembly would meet the definition that the cells within the assembly are firing at a high rate in the same time frame. Successions of assembly activation could be the basis of how to form thoughts and give a mental representation for different concepts. This idea was introduced by Donald Hebb in The Organization of Behavior. The assembly will eventually shut down if the source of stimulation ceases or a competing assembly can provide inhibition. This theory can be extended to explain how learning occurs and possible phenomena such as perspective-rivalry. For instance, in optical illusions, such as the nectar cube shown below, there are two different possible interpretations of the perspective. It is possible that two different assemblies represent the different perspectives. In order to switch perspectives one must shut down the current perspective assembly and activate the alternate perspective. There could also be a hierarchy of assemblies. Sub-assemblies could give rise to larger assemblies. Due to synaptic plasticity different sub-assemblies that are activated in similar orders through repetition allows new larger Assemblies to form. If sufficient clues that are represented by smaller sub-assemblies are activated, the brain will be able to activate the larger assembly.

The Neckar Cube has two possible perspectives. Cell Assembly rivalry could be responsible for competing interpretations.

Module Goal To get some insight on whether this is plausible we want to construct biologically accurate computational models and test assembly activation and inhibition. This module will follow the paper Modeling Hebbian cell assemblies comprised of cortical neurons, by A Lanser and E Fransen. It will show how to reproduce the results contained in the paper and the mathematics involved. Here is a brief summary of what is ahead: 1.) Go over the neuron models and the different ion channels. 2.) Show how synaptic input is modeled. 3.) Show how differential equations are solved with time stepping. 4.) Show how patterns which represent different assemblies are used to "train" a network and produce weights. 5.) Give parameters that are used and explain parameters that are dependent on simulations. 6.) Run simulations by going through the matlab prompt and display results.

Cell Models Below is a quick description of the basic dynamics of neurons and how messages are relayed in neuronal networks. The mathematical model that is used in the simulations is given in detail. The basic function of a neuron is the ability to relay information by means of an action potential. This is a quick depolarization burst which originates at the cell soma (main body) and travels down the axon reaching boutons. At a bouton the depolarization will cause neurotransmitters to be released into the synaptic cleft, a small space between the bouton and a spine of the target cell. The neurotransmitter then opens ion channels in the postsynaptic receptor area that can serve to excite or depress the post-synaptic neuron, depending on whether the sending neuron is excitatory or inhibitory. A cell achieves an action potential when it receive enough excitation from other cells. When the cell voltage reaches a certain threshold its firing mechanisms kick in and generate an action potential. The firing mechanism consists of voltage gated ion channels in the cell's soma membrane. These channels permability are dependent on voltage levels and are specialized to only allow specific ions to pass through. The standard action potential starts with stimuli which raises voltage sufficiently to open Na+ channels. Then Na+ floods in (due to diffusion) causing a depolarization of the cell. When the voltage reaches the near peak depolarization K+ gating channels open allowing K+ to rush out, lowering the potential back to rest. There are additional ions which have a role in cell dynamics. One important ion that helps regulate firing rates is Ca2+. Ca2+ levels are important because they open Ca2+ dependent K+ channels. More K+ channels lead to greater repolarization, increasing refractory times between firing and leading to lower activity. We will look at two different sources of Ca2+. One, through voltage gated channels. The second from NMDA receptor channels, which will be discussed in further detail ahead.

Cell Types There are two basic types of cells: excitatory and inhibitory. When excitatory cells fire they increase the potential in the neurons they synapse onto. Whereas, inhibitory decrease the potential. For our purposes we will use different parameters to the two cells different characteristics.

Mathematical Model The excitatory cell model used has 4 compartments per cell. The first compartment represents the soma. The other three are a chain of dendrites, which represent a collapsed dendritic tree at different distances from the soma. In the case of an inhibitory cell there will only be 2 total compartments. Cells will synapse onto compartment 4 if it is a E (excitatory) to E connection. Compartment 1 if it is an I (inhibitory) to E connection. Compartment 2 if it is an E to I connection. Below is the main differential equation that governs the cell voltage in a compartment. The first term in the fraction represents the leak current. The second term is the summation of current that comes in from neighboring compartments. The third term represents Ion channel currents; they are only present at the cell soma. The fourth term is current entering from synaptic input. d V d t = V l e a k - V G m + ∑ V c o m p - V G c o r e + I c h a n n e l s + I s y n C m Variable Description V Compartment Potential E l e a k Outside Potential G m Membrane conductance E c o m p Neighbor compartment voltage G c o r e Conductance between Compartments I c h a n n e l s Ion channel current (only in soma compartment) C m Membrane capacitance

Ion Channels The cell dynamics use a model that is derived from the Hodgkin-Huxley model. The ion currents take the general form of a driving force (which is due to the different concentrations of the ions on the inside and outside of the cell) being multiplied by a conductance term (the reciprocal of resistance) and a product of voltage dependant variables that indicate the fraction of the ion channels open. For example, the sodium current is: I N a = ( V N a - V s o m a ) G N a m 3 h . Where VNa is the reversal potential (indicating the concentration of Na+ in the surrounding fluid). GNa is the conductance of Na+. m is called the activation variable of the Na+ channels and h is the inactivation variable. They are governed by the following differential equations: d m d t = α m ( 1 - m ) - β m m Here αm is the opening rate and βm is the closing rate. These rates are dependent on the soma potential. A, B and C are fitting parameters that will later be prescribed. α m = A ( V s o m a - B ) 1 - e ( B - V s o m a ) / C β m = A ( B - V s o m a ) 1 - e ( V s o m a - B ) / C The inactivation variable h is slightly different, and has its own fitting parameters. d h d t = α h ( 1 - h ) - β h h α h = A ( B - V s o m a ) 1 - e ( V s o m a - B ) / C β h = A 1 + e ( B - V s o m a ) / C The K+ channels are similar. There is only one gating variable n. I K = ( V K - V s o m a ) G K n 4 d h d t = α h ( 1 - h ) - β h h α h = A ( V s o m a - B ) 1 - e ( B - V s o m a ) / C β h = A ( B - V s o m a ) 1 - e ( V s o m a - B ) / C For more details on the Hodgkin Huxley model visit http://www.caam.rice.edu/caam415/ or look at the series of papers this module is designed after. The change in net potential due the calcium flow is negligible. However, we still need to track it in order to control the Calcium dependent Potassium channels. As cell repeatedly fires calcium will build up leading to stronger repolarizations this will help allow for firing rate adaptation and provide the assembly a mechanism through which it will die out eventually. I q = ( V q - V s o m a ) G q n 4 d q d t = α q ( 1 - h ) - β q q α q = A ( V s o m a - B ) 1 - e ( B - V s o m a ) / C β q = A ( B - V s o m a ) 1 - e ( V s o m a - B ) / C The Calcium dependent Potassium channels: I K ( C a ) = ( V K - V S o m a ) G K ( C a ) [ C a A P ] d [ C a A P ] d t = ( V C a - V s o m a ) ρ A P q 5 - δ A P [ C a A P ] [CaAP] is the intracellular calcium level. Constants ρAP and δAP are the rates of calcium ion influx and efflux (decay). The second form of Calcium comes through NMDA channels that are located on the post-synapse. It is important to separate this Calcium source since it enters and leaves the cell with different time constants due to the different method it enters the cell. This contribution will be discussed in more details in the next section. The Calcium pool [CaNMDA] will be added to the total Ca2+ count changing the total current from Calcium dependent Potassium channels to: I K ( C a ) = ( V K - V S o m a ) G K ( C a ) ( [ C a A P ] + [ C a N M D A ] )

Synaptic Connections Cells communicate with synaptic connections where the bouton and post-synapse meet. The firing cell releases neurotransmitter into the synaptic cleft. These chemical open up ion channels in the postsynaptic receptor cell. This allows either a gain in postsynaptic potential (if the presynaptic cell is excitatory) or a depression in potential (for a inhibitory cell). The model for the synaptic current is show below: I s y n = ( V s y n - V ) G s y n s . Here s is an activation level variable. It will only take value 0 or 1. After a presynaptic cell fires the s variable will be held at 1 for fixed number of milliseconds, then reset to 0. The excitatory receptor area will consist of AMPA receptor and NMDA receptors. The AMPA receptor takes the form above and Vsyn would typically be around 0 mV for an excitatory connection. An inhibitory connection will have a reversal potential, Vsyn, below or near rest (around -85 mV). The second type of receptor Channel is the NMDA channel. These receptor channels are blocked by Mg2+ when the post-synaptic cell in near resting potential. Only come out when the receptive cell has been depolarized will the Mg2+ pop out allowing Na+ and Ca+ to pass through. This allows so that cells can continue firing once they have been activated.These channels are dependent on neurotransmitter from the sending cell and the post-synaptic cell firing. This feature helps for LTP (long term potentiation). The Mg2+ is modeled by a voltage dependant ODE between 0 and 1, p. The form for this current and it's gating variable p is the following: I N M D A = ( V N M D A - V ) G N M D A p s . d p d t = α p ( 1 - p ) - β p p . α p = A e E / C β p = A e - E / C The NMDA channels are also permeable to Ca2+ so we need to track another pool of calcium. This calcium is modeled by a differential equation. The influx rate is dependent on the s and p variable, as well as a characteristic influx conductance term, ρNMDA. This influx term varies from each synapse and is proportional to the synapse strength (something which is determined in later a section when connection weighting is addressed). The efflux is a constant δNMDA. The s variable here is the sum of all the binary s variables which connection to the specific cells whose Calcium pool is being calculated. d [ C a N M D A ] d t = ( V N M D A - V ) ρ N M D A p s - δ N M D A [ C a N M D A ] This term actually looks a little more complicated. There are generally many synapses where Ca2+ can enter a neuron and the entry rate at each synapse is a product of the activation variable skj (which is 1 if the pre-synaptic neuron j which synapses onto k fired recently, 0 if not) and the connection weight ρNMDA, which is proportional to the conductance wieghting (i.e. ρNMDA=NMDAscalarGkj). We must write this term as a summation of all the Ca2+ coming in from different synapses. There are now subscripts and superscripts j and k being used to indicate which cells are be tracked. In the below equation is the differential equation for cell k. The j index is used for other cells that synapse onto k. The compartment here is assumed to be the soma as we will not bother to track the calcium flow from the end compartment to the soma. d [ C a N M D A ] k d t = ( V N M D A - V k ) p k ∑ j onto k N M D A s c a l a r G k j s k j - δ N M D A [ C a N M D A ] k

Computational Methods for Solving Differential Equations There are several differential equations that need to be solved in order to advance the system. A time-stepping method is used and updates each variable given the selected step size (typically .01 milliseconds). A hybrid Euler scheme is used. In a differential equation the derivative is replaced with the approximation (vi+1-vi)/dt. The right hand side vi+1 is plugged in for v when the equation can be explicitly solved for, this is backward Euler which has better stability properties. However, when we cannot solve for vi+1 we use vi instead, this is a forward Euler stepping scheme. In our main differential equation for cell voltage some terms we can use the backward method. The m,h,n,q,[CaAP] and [CaNMDA] variables have to updated using a forward Euler update first. Then they are plugging into their ion current components of the main cell voltage differential equation. This is all then solved for Vi+1. The set up is below. Remember the ion currents are only at the soma, below we imagine we are updating the soma compartment. Other compartments are done the same way except the ion currents are not there and some parameters such as Cm take different values. We start with an example of how the m gating variable would be updated. ( m i + 1 - m i ) d t = α m [ V i ] ( 1 - m i + 1 ) - β m [ V i ] m i + 1 = A α m ( V i - B α m ) 1 - e ( B α m - V i ) / C α m ( 1 - m i + 1 ) - A β m ( B β m - V i ) 1 - e ( V i - B β m ) / C β m m i + 1 Solving for mi+1 gives: m i + 1 = m i + d t α m [ V i ] 1 + d t α m [ V i ] + β m [ V i ] Update n,h,q,[CaAP],[CaNMDA] likewise solving for the i+1 step. Notice αm and βm are functions of voltage. Plug these results into below equation which is then solved for Vi+1 V i + 1 - V i d t = V l e a k - V i + 1 G m + ∑ V c o m p - V i + 1 G c o r e + I N a + I K + I C a + I K ( C a ) + I s y n C m The ion channels have the Vi+1 term in them. For example the Na+ channel looks like this: I N a = ( V N a - V i + 1 ) G N a m i + 1 3 h i + 1 The final result term after solving for Vi+1: V i + 1 = C m V i + d t ( G m V l e a k + G N a m i + 1 3 h i + 1 V N a + G K n i + 1 4 V K + G K ( C a ) V K ( [ C a A P ] i + 1 + [ C a N M D A ] i + 1 ) + I S y n + I S t i m C m + d t ( G m + G N a m i + 1 3 h i + 1 + G K n i + 1 4 + G K ( C a ) ( [ C a A P ] i + 1 + [ C a N M D A ] i + 1 )

Building Networks After constructing models for the individual neurons it is necessary to connect the cells. The mission is to have a small network (50 excitatory and 50 inhibitory neurons will be used in the simulations, more details on why in a bit) and have the cells trained with pre-selected assemblies/patterns. Given different patterns or events we want to construct a weighting scheme that will assign a proper weighting. Proper in the sense that cells within a pattern will have strong excitatory connections and cells that are ever in the same pattern will have inhibitory connections. The goal is to have the network able to activate patterns. This means that when a sufficient number of cells belonging a pattern are stimulated it will fully resolve. The weighting algorithm will output weights that are effective in proportion to another, but will require proper scaling. Should the overall weighting magnitudes be too high the network will seizure (meaning all the cells start firing). If the weighting magnitudes are too low, the network will go dormant once stimulation is removed.

Picking Weights The term weight is used loosely; it signifies the choice of the measure of conductance at synapse where two neurons meet. This parameter will suffice as way of quantifying the strength of a synapse which is physiologically determined by the amount of neurotransmitter that would be released into the cleft and the number of receptor channels. Look back at the form for a synaptic current to see the conductance term we now will give more detail on: I s y n = ( V s y n - V ) G s y n s . The Gsyn is the conductance term that reflect the weighting. It is expressed as the product of a weighting term that is an output of the weighting algorithm and an overall weighting scale constant that will determined later. Thus, Gij=wij*weightsyn. Note: there will be several different weightsyn constants for the several synaptic currents used in the model. The scheme used to produce weights will use a probabilistic approach. Using a simple intuitive approach we want to characterize the level of support a cell gets by the term: s q = b q + ∑ h w h q π h Here sq is the support level of the cell. The term bq will represent a bias of a cell (this would be a way of characterizing a cells firing property independent from other cell activity, since some cells may appear simply more activate in general). whq is the connection weight from h to q. Variable πh is the activity level of cell h. Assume either 0 or 1 then we reduce this system to a summation over the active cells in all the different patterns, set A. s q = b q + ∑ h ϵ A w h q From conditional probability rule p ( q | i ) = p ( q & i ) p ( i ) arises Bayes Rule, which can then be iterated. p ( q | i ) = p ( q ) p ( q | i ) p ( i ) . p ( q | i & j & k & . . . ) = p ( q ) p ( q | i ) p ( i ) p ( q | j ) p ( j ) p ( q | k ) p ( k ) . . . = p ( q ) ∏ h ϵ A p ( h | q ) p ( h ) . Taking the logrithm of this term gives us something in the form we wanted: l o g p ( q | A ) = l o g p ( q ) + ∑ h ϵ A l o g p ( h | q ) p ( q ) . Thus we have that: w h q = l o g p ( h | q ) p ( q ) = l o g p ( h & q ) p ( q ) p ( h ) . The method for computing the probabilities is now explicated. For single probabilities (such as p(q)), we simply look at each pattern and compute the proportion in which cell q is are active. If either p(q) or p(h) be zero set whq to zero. For joint probabilities (such as p(h&q)), compute the proportion in which both are active. If both p(q) and p(h) are not zero but p(h&q) is zero, set whq to log(1/#patterns). There exists a more complicated version in which different events/patterns carry more significance than others and are then weighted accordingly by κ(α), this capability is not implemented in our simulations but the code can handle this option. The code for building the weight matrix (where entry wij is connection strength from i to j) is carried out by the code make_W.m. The method that produces the weight matrix W will produce a symmetric matrix. This should not be that limiting on the network model, but is not exactly representative of physiology.

Hooking up the Network The output from the weight matrix will yield positive and negative entries. If the entry is positive we have an excitatory connection. If an entry is negative there is an inhibitory connection. If the absolute value is less than a chosen tolerance then there is no connection. Now comes an important step that is not entirely biologically correct but serves to make the over network model simpler. The patterns that we used in constructing the weight matrix in fact will consist of all excitatory cells. For each of these excitatory cells we will now designate a inhibitory companion cell through that only synapses onto the excitatory cell. All inhibition of the excitatory cells is achieved by other excitatory cells that synapse onto the inhibitory companion. Thus, there are an equal number of E and I cells. (as mentioned earlier there will be 50 E and 50 I cells in the simulations produced ahead).

Here we see an 8 cells network that has a the same conventions. The excitatory cells have 4 compartments each. Excitatory cell connections synapse onto the end compartment of other E-Cells and the second compartment of I-Cells (Since I-Cells only have 2 compartments). Inhibitory cells synapse onto the Soma's of the excitatory cells.

This is a diagram for an excitatory-inhibitory pair. In some sense we can regard this as a full neuron unit. Here the compartments are not show as they are above.

Simulations The bulk of the theory is completed. Now, the specifications for the simulations are presented. The parameters that are used in the biological model presented earlier are given values. Then there is a few examples of how to use the matlab code. The code used can be downloaded at http://www.caam.rice.edu/~cox/hippo.html#code .

Network Size and Patterns In the simulations there will be 50 E-cells and 50 I-cells. The Network will be "trained" to encode assemblies by running the weight producing algorithm on a set of patterns, each which represents a cell assembly (remember weights are only produced between the E-cells initially). There will be 8 patterns that contain 8 E-Cells in each. Each "event" or pattern will carry equal significance when the weighting algorithm runs. To help visualize the patterns look at the diagram below.

Parameters For Cells Below are tables with parameter that are given exactly from A Lansner and E Fransen. Parameter E-Cell I-Cell Vleak (mV) -50 -70 G c o r e 0.04 0.0638 Gm Soma (μS) 0.0032 0.0016 Cm Soma (nF) 0.032 0.016 Gm(μS) Dentrites 0.0096 0.0096 Cm Dentrites (nF) 0.288 0.288 V N a 40 50 G N a ( μ S ) 1.0 1.0 V K -70 -90 G K ( μ S ) 0.5 1.0 V C a 150 150 G C a <