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.

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:

Where VNaVNa is the reversal potential (indicating the concentration of Na+Na+ in the surrounding fluid). GNaGNa is the conductance of Na+Na+. mm is called the activation variable of the Na+Na+ channels and hh is the inactivation variable. They are governed by the following differential equations:

Here αmαm is the opening rate and βmβ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.

The inactivation variable hh is slightly different, and has its own fitting parameters.

The K+K+ channels are similar. There is only one gating variable nn.

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.

The Calcium dependent Potassium channels:

[CaAP][CaAP] is the intracellular calcium level. Constants ρAPρAP and δAPδ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][CaNMDA] will be added to the total Ca2+Ca2+ count changing the total current from Calcium dependent Potassium channels to:

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:

Here ss is an activation level variable. It will only take value 0 or 1. After a presynaptic cell fires the ss 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 VsynVsyn would typically be around 0 mV for an excitatory connection. An inhibitory connection will have a reversal potential, VsynVsyn, below or near rest (around -85 mV).

The second type of receptor Channel is the NMDA channel. These receptor channels are blocked by Mg2+Mg2+ when the post-synaptic cell in near resting potential. Only come out when the receptive cell has been depolarized will the Mg2+Mg2+ pop out allowing Na+Na+ and Ca+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+Mg2+ is modeled by a voltage dependant ODE between 0 and 1, pp. The form for this current and it's gating variable pp is the following:

The NMDA channels are also permeable to Ca2+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 ss and pp variable, as well as a characteristic influx conductance term, ρNMDAρ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δNMDA. The ss variable here is the sum of all the binary s variables which connection to the specific cells whose Calcium pool is being calculated.

This term actually looks a little more complicated. There are generally many synapses where Ca2+Ca2+ can enter a neuron and the entry rate at each synapse is a product of the activation variable skjskj (which is 1 if the pre-synaptic neuron j which synapses onto k fired recently, 0 if not) and the connection weight ρNMDAρNMDA, which is proportional to the conductance wieghting (i.e. ρNMDA=NMDAscalarGkjρNMDA=NMDAscalarGkj). We must write this term as a summation of all the Ca2+Ca2+ coming in from different synapses. There are now subscripts and superscripts jj and kk being used to indicate which cells are be tracked. In the below equation is the differential equation for cell kk. The jj index is used for other cells that synapse onto kk. 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.

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(vi+1-vi)/dt. The right hand side vi+1vi+1 is plugged in for vv when the equation can be explicitly solved for, this is backward Euler which has better stability properties. However, when we cannot solve for vi+1vi+1 we use vivi 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]m,h,n,q,[CaAP] and [CaNMDA][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+1Vi+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 CmCm take different values. We start with an example of how the mm gating variable would be updated.

Solving for mi+1mi+1 gives:

Update n,h,q,[CaAP],[CaNMDA]n,h,q,[CaAP],[CaNMDA] likewise solving for the i+1i+1 step. Notice αmαm and βmβm are functions of voltage. Plug these results into below equation which is then solved for Vi+1Vi+1

The ion channels have the Vi+1Vi+1 term in them. For example the Na+Na+ channel looks like this:

The final result term after solving for Vi+1Vi+1:

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:

The GsynGsyn 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*weightsynGij=wij*weightsyn. Note: there will be several different weightsynweightsyn 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:

Here sqsq is the support level of the cell. The term bqbq will represent a bias of a cell (this would be a way of characterizing a cell’s firing property independent from other cell activity, since some cells may appear simply more activate in general). whqwhq is the connection weight from hh to qq. Variable πhπh is the activity level of cell hh. Assume either 0 or 1 then we reduce this system to a summation over the active cells in all the different patterns, set A.

From conditional probability rule

arises Bayes Rule, which can then be iterated.

Taking the logrithm of this term gives us something in the form we wanted:

Thus we have that:

The method for computing the probabilities is now explicated. For single probabilities (such as p(q)p(q)), we simply look at each pattern and compute the proportion in which cell q is are active. If either p(q)p(q) or p(h)p(h) be zero set whqwhq to zero. For joint probabilities (such as p(h&q)p(h&q)), compute the proportion in which both are active. If both p(q)p(q) and p(h)p(h) are not zero but p(h&q)p(h&q) is zero, set whqwhq to log(1/#patterns)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 wijwij is connection strength from ii to jj) is carried out by the code make__W.m.

The method that produces the weight matrix WW 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).

Figure 2: 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.

Figure 3: 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.

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.

Remaining Parameters

There are a few parameters that remain to be given values. A Lanser and Fransen do not give exact values for these unknown parameters, but instead a desired range for EPSP's (Excitatory Post Synaptic Potentials) that the parameters help determine. By experimenting with a single neuron we can find the ideal range for these unknown parameters. Hence, we can then scale the resulting Weights from the algorithm so they are mapped to the ideal range. We need an overall weighting scale constant for connections from E to E cells (AMPA), wE2EwE2E. A weighting scale constant for E to I cells, wE2IwE2I Also, a scale general conductance constant for all I to E connections, gI2EgI2E. The NMDA connections will be proportional to the AMPA weight but scaled by a constant wNMDAscalarwNMDAscalar. The influx parameter ρNMDAρNMDA for the [CaNMDA][CaNMDA] pool is also proportional to the AMPA weight. Lastly, the binary variable ss needs to be given a time duration to stay active for. This time duration may be chosen different for each type of connection class as well (meaning E to E, E to I, and I to E).

Most of these unfixed variables relate to the strength of connections between cells and the EPSP’s or IPSP’s they create. The reason they won't be fixed for all different size simulations is that the values should vary according to size of the networks and assemblies. For instance, if we model on large networks that have large assemblies we would desire that a cell would require synaptic input from more cells in order to fire, than when the network is small. By picking these numbers accordingly we will be able to roughly determine how many cells of an assembly need to be firing in order to fully activate the assembly.

Simulation 1

I have created a code called pattern__maker.m which is an easy way to get a matrix will a designated number of cells, number of patterns, and number of cells per pattern. We will start with a network with 50 excitatory cells and train it with 8 patterns, containing 8 cells in each pattern. Type help pattern__maker into a matlab prompt to get the details:

>> help pattern_maker
pattern_maker.m - creates a pattern matrix. Each row represents a
pattern with 1's for active cells, 0's for inactive cells.
pattern = pattern_maker(n_cells,n_patterns,n_active)
where: n_cells = number of cells
n_pattern = number of patterns
n_active = number active units per pattern
output: pattern = n_patterns x n_cells matrix with each row a
randomized pattern with n_active units

Now lets create a variable pp which will store the pattern matrix.

>> pp = pattern_maker(50,8,8);

Now that we have a pattern we can run the main piece of code that will do the simulation: LF__network.m. Note that many of the parameters are only accessible in the code and would have to be directly edited if the user wishes to change them.

>> help LF_network
LF_network.m - This function runs the Lansner-Fransen type simulations as
detailed in their paper - "Modelling Hebbian cell assemblies comprised of
cortical neurons". There are an equal number of excitory and inhibitory
cells. Each excitory cell has a companion inhibitory cell which synapses
onto it. The network is trained with different patterns loaded in. Then
its calls a code which produces connection strengths. Positive strengths
form Excitory to Excitory connections. Negative strengths form
connections from Excitory to Inhibitory cells.
[v_ex,v_in] = LF_network(dt,Tcutoff,Tfin,I0,stim_cells,Patterns)
where: dt = time step
Tcutoff = time cutoff for stimulated cells
Tfin = time end simulation
I0 = stimulation size for cells
stim_cells = vector of cells recieving stimulus
Patterns = matrix where each row is a pattern the network is
trained with
returns: v_ex = Voltages for excitory cells, indexed as
ve(cell,timestep,compartment)
v_in = Voltages for inhibitory cells, indexed as
vi(cell,timestep,compartment)
ex. pp = pattern_maker(50,8,8)
[v_ex, v_in] = LF_network(.01,50,350,1,[1 6 8 10 25],pp);

As explained in the help file above one must choose the cell numbers to stimulate. In order to pick several cells in one assembly we will enter the following command to get the indices of one of the patterns:

>> find(p(1,:)==1)
ans =
19 23 28 29 32 35 36 42

Thus we will choose to stimulate cells 19, 23, 28, 29 and add the noise stimuli of cells 45, 46, and 47. The simulation will last for 350 ms and the stimulation will last for 50 ms. A 1.5 nA stimulus is injected into each respective cell soma. After the output is saved we load the results in a code which displays the results: viewer.m.

>> [v_ex,v_in] = LF_network(.01,50,350,1.5,[19 23 28 29 45 46 47],pp);
>> help viewer
viewer.m takes in output from LF_network and plots voltage traces of all
the cells. If you wish to only look at one cell enter that as a
optional 3rd input.
viewer(v_ex,v_in,cell_trace)
>> viewer(v_ex,v_in)

Figure 4: Cells 19, 23, 28 and 29 of pattern [19 23 28 29 32 35 36 42] are stimulated as well as noise cells 45, 46, and 47. The blue traces are the excitatory cells, red are the inhibitory companions. The full assembly activates and then dies out around 250ms. The noise cells shut down before 50ms even passes.

We can get a better look by looking up close at the trace of cell 36, one of the cells in the assembly that did not receive a stimuli. Just type viewer(v__ex,v__in,36) into the prompt.

Figure 5: This is a close up of cell 36 Excitatory and Inhibitory from the simulation above. The blue trace of the excitatory cell slowly increases until it reaches a threshold where upon it fires rapidly. The interspike intervals decrease until the cell shut off. The inhibitory cell show by the red trace receives no input once the noise cells have been shut down.

Figure 6: Here the two Ca2+Ca2+ pools are plotted for excitatory cell 36. Note the AP calcium enters quicker, the NMDA calcium has a slower time course and helps shut down the assembly.

Simulation 2 - Pattern Competition

In this simulation we will stimulate cells in two different patterns and see if one will "win" by activating it’s full pattern and suppressing the competing pattern. We will use the pattern from above and also another pattern.

>> find(p(6,:)==1)
ans =
11 12 19 22 27 29 33 39

This pattern 6 should be interesting since it has two overlapping cells with pattern 1 (that is ## 19 and ## 29). Let’s run a simulation that picks one of the overlapping cells (19), 4 others from pattern 1 (28 32 35 36) and 3 more from pattern ## 6 (11 12 22).

>> [v_ex,v_in] = LF_network(.01,50,350,1.5,[11 12 19 28 32 35 22 36],pp);

Figure 7: 5 cells from pattern 1 and 4 cells from pattern 6 are stimulated. Pattern 1 resolves and shuts down around 250 ms

We can now try stimulated 4 cells from pattern #1#1 and 5 cells from pattern #6#6.

>> [v_ex,v_in] = LF_network(.01,50,350,1.5,[11 12 19 22 27 28 32 35 ],pp);

Figure 8: 4 cells from pattern 1 and 5 cells from pattern 6 are stimulated. This time Pattern 6 resolves and shuts down around 250 ms