Abstract

I. Introduction

Developing better methods for determining depth of anesthesia (DOA) requires good understanding about brain functions during general anesthesia. Mean filed models (MFM) because of their remarkable properties of handling thousands of similar neurons in one population are one of the best candidates to study brain functions and activities during anesthesia [1], [2].

Recently we have proposed an enhanced version of Bojak MFM by adding a slow modulatory mechanism of neural firing rate in it [3]. This mechanism was added to Bojak model based on some real neuro-physiological experiments [4],[5]. Here, we show how the role of fast and slow variables in our enhanced MFM (EMFM) can be investigated. To do this, first of all we briefly introduce our EMFM and then adiabatically simplify the model to have only two fast (inhibitory and excitatory membrane potentials) and one slow variables. Then we obtain isoclines of the slow variable and excitatory membrane potential that is assumed to be proportional to EEG signal. To be able to investigate on the role of slow and fast variables from waking to deep anesthesia, we sketch the vector field of these variables. Since we have two fast variables and one slow variable, the vector field would be three-dimensional. Here, we have proposed a method that combines the two fast variables into one variable. This enables us to study the role of variables of the model by a two-dimensional vector field that is sketched on the isoclines plane.

This method was used to investigate the role of slow modulatory variable in our EMFM when a typical anesthetic drug concentration is changed from zero to a certain amount that produces a very deep anesthesia. Investigations show that the slow variable can generate a bi-stable neural firing. This bi-stability that occurs in delta frequency range (1-4Hz) is the main reason that brain delta rhythms can be recorded on scalp with high amplitudes much greater than those in alpha or beta rhythms during sleep or deep anesthesia.

II. Basics of our enhanced mean-field model

Our EMFM has been designed based on Bojak model [2]. Figure 1 illustrates a schematic diagram of the basic and enhanced MFM. In EMFM, brain is considered as a homogeneous media, so its activities are studied by identical macrocolumns. Macrocolumns are consisted of two neuronal populations, excitatory and inhibitory. Each macrocolumn is represented by 9 connected differential equations and it would be the whole differential equations of the model in homogeneous case. These equations have been brought briefly below (see [3] and [2] for more details):

Open in a separate window

Figure 1

Schematic diagram of the inter/intra connections in the macrocolumns. The switch determines which model has been selected, ‘basic’ or ‘enhanced’.

Equation (1) expresses mean membrane potentials of inhibitory and excitatory populations by flowing currents to and from an RC circuit. Voltage sources of this RC circuit are: leakage ( $h_k^{\mathit{rest}}-h_k$), excitatory synapses Ψ_ek(h_k)I_ek(h_e) and inhibitory synapses Ψ_ik(h_k)I_ik(h_i). I_jk is the voltage response of population j to k -type spike rate when it is in resting potential. Ψ_jk(h_k) is used as a compensation factor for resting potential. Equation (2) describes post synaptic potential of population j by G_j, γ_j and _j in a second order differential equation, k -type spike rate is composed of three kinds of sources: locally in a same macrocolumn (S_j), distant from other macrocolumns(Φ_jk) and subcortical input noises (p_jk). $N_{jk}^{\beta }$ is the number of connections from population j to population k within a macrocolumn and $N_{jk}^{\alpha }$ is the same quantity but it is used for inter connections of macrocolumns.

In EMFM, there is also a slow mechanism that modulates the excitatory firing rate. This slow mechanism is used as a simple representation of many kinds of slow intrinsic ionic channels on neural membranes that have not been considered directly in the RC circuit of neural cells.

We have assumed that a slow process that is called s modulates the firing rate of excitatory population, s follows s_∞ based on the first order differential equation in (5). τ_s is a high value time constant and s_∞ is a descending sigmoid function of the excitatory membrane potential and represents the tendency of a cell to be in up state. It is meaningful when the potential is decreased there is more tendency for neural cells to be in up state.

In the basic MFM, excitatory firing rate S_e is mainly represented by excitatory mean potential h_e (see Eq. (7)) but in EMFM, s also influences the value of S_e.

Equation (6) has gathered two different terms of firing rates, S_e(h_e) and $S_e^{max}$. Normalized weighted of these two terms determines the overall firing rate of excitatory populations. Weighting coefficients F₁(s) and F₂(s) are anti-symmetric sigmoid functions of s. B which lays between zero and one is a free parameter that determines the weighting balance of each firing term. θ_F and g_F respectively determine the turning point and slope of sigmoid functions. g_F must take a value that is smaller than zero but θ_F takes a value between zero and one.

Administration of anesthetic drugs varies the shape of PSP functions. The effect of anesthetic drug concentration can be applied on EMFM by declaring γ_j, _j and G_j as functions of drug concentration. These parameters mainly determine PSP by equation (2). Generally a Hill equation is used to express this parameter as functions of drug concentration [2].

III. Model solutions

B. Isoclines of slow mechanism and excitatory membrane potential

One of the best ways to study the influence of slow variable s on the network in different anesthetic concentration is to sketch ds/dt = 0 and dh_e/dt = 0 isoclines in a same plane (s h_e isoclines) [6]. In the case of stable equilibrium state, s is reached its target value s_∞ but in an unstable mode s and s_∞ may have different values due to the high value of τ_s. This can produce some special network activities that can be investigated and interpreted by isoclines.

To sketch s h_e isoclines any variables in the nine connected differential equations is set to its equilibrium value. Suppose that from (5) a set of ŝ values that correspond to a given set of stationary h_e values, have been calculated and plotted versus h_e. For each point on the obtained curve we have ds/dt = 0, so this curve is called ds/dt = 0 isocline (simply isocline #1). Reversely it is possible to obtain ĥ_e as a function of equilibrium s values to obtain dh_e/dt = 0 isocline (isocline #2). To do this, s is altered between zero to one (its valid range) with a desired step size (e.g. 0.05) and corresponding with each s value in this range, equilibrium values of h_e are obtained using the root search method described previously. Now, if we properly superimpose these two isoclines as it has shown in Figure 2, intersections of them indicates the equilibrium solutions of the differential equations that were found previously because ds/dt = dh_e/dt = 0 for the intersected points. Any other point on isoclines plane may have positive or negative values of ds/dt and dh_e/dt. These values are represented by a vector field on the isoclines plane. Vectors lengths and directions show how and to what extend, s and h_e are changed if the network state is positioned on them, so they can predict the trajectories. It is assumed that comparing to s, h_e and h_i, other variables in the model are enough fast and reach their steady state rapidly. This assumption reduces our model to a three-dimensional model; as a result we can trace the trajectories by a 3D vector field.

Open in a separate window

Figure 2

A typical dh_e/dt = 0 and ds/dt = 0 isoclines and six candidate planar vectors that are originated from ( $s=0.3,h_e^{•}\pm 1mV$)

IV. Vector field and network behavior

In this section we show how we can obtain a two-dimensional vector field on the isoclines plane that can predict the behavior of the numerically simulated signals (trajectories). For this purpose, we illustrate two s h_e isoclines plane corresponding with two different drug concentrations that exhibit different dynamics. Using this method, we will able to see how the slow mechanism can aperiodically switch neural state between up and down when drug concentration is increased.

A. Obtaining the 2D vector field on isoclines plane

Figure 3A illustrates s h_e isoclines (dashed and dotted lines) and a 5sec superimposed trace of s(t) and h_e(t) signals when drug concentration is set to zero. Time course of h_e(t) signal has been illustrated in Figure 3B. Arrows in Figure 3A represent the planar vector field for some parts of isoclines plane. For each s between 0 and 1 in 0.05 steps, these vectors have been calculated in $h_e^{•}+\mathrm{\Delta }{h}_e$ where Δh_e = {−1,1} mV and $h_e^{•}$ are red dots on isocline #2. Calculating vertical component of these vectors (in s direction) requires inserting $h_e^{•}+\mathrm{\Delta }{h}_e$ in (5) and obtaining ds/dt from this equation. Horizontal components are obtained from (1) where k ≡ e. But since in this equation, h_i also influences the value of dh_e/dt, to be able to study the response of the model by vector field, we have to obtain a 3D vector field (dh_e/dt, dh_i/dt, ds/dt) that makes the problem a little more complicated. Instead of using a 3D vector field, we have proposed a method that reduces the vector dimensions to two. To do this, for each point located in ( $s,h_e^{•}+\mathrm{\Delta }{h}_e$) on isoclines plane, we tested some $h_i^{•}+\mathrm{\Delta }{h}_i$ values where Δh_i was changed from −1mv to +1mv. Our experiments show that a 2D vector field (dh_e/dt, ds/dt) properly traces the simulated signals in the isoclines plane if Δh_i κ Δh_e where $\kappa =h_i^{•}/h_e^{•}$. Although using this method vectors are coarsely obtained, it is possible to match all traces of simulated signals with vector field directions. This approximation can be used as a useful method for tuning the model parameters and investigating on the role of a specific variable in the model.

Open in a separate window

Figure 3

s h_e isoclines and their corresponding vector field and five seconds of s(t) and h_e(t) trajectories for zero drug concentration (A). The time course of h_e(t) trajectory (B).

B. Studying the role of the slow mechanism in EMFM Using the 2D vector field

In all of our simulations we used Euler method with 0.1ms step size and a down sampler to generate 10-second signals having 400Hz sampling frequency.

Here we show the planar vector field is able to predict the model behavior before performing any numerical simulation. It can be inferred from isocline #2 of Figure 3A that regardless the value of s there is typically one equilibrium solution about high h_e potentials that expresses being in the up state. Even for s’s that have three different equilibrium points on isocline #2 (0 < s < 0.4 ), the highest potential equilibrium point (up state) has the most occurrence probability. This is evident from the vector near the three equilibrium points. Upward short length vectors about the lowest potential equilibrium point (down state) indicate that there is not enough attraction to this point and the simulated trace cannot reach to this point because there is much distance between the equilibrium point in up and down states. h_e cannot also be settled in the middle equilibrium point because the vectors about this point go away from it (unstable point). As a result, h_e and other model parameters are settled in their up state position.

In Figure 4A the network balance has been changed toward more inhibition by increasing anesthetic drug concentration. It can be inferred from this figure that the right branch of isocline #2 has moved a little to more negative potentials also the U-turn has moved to higher s values. Note that s positively modulates the excitatory firing rate in (6); so moving the U-turn in positive s direction indicates that the balance has moved to inhibition because more extra modulation of excitatory firing rate is required to have only one equilibrium point on isocline #2 (this happens when s ≈ 0.55 ). It means for an identical s value that isocline #2 of zero drug concentration just has one equilibrium point; here it may have three equilibrium points. One of the effects of shifting the right branch of the isocline #2 and reducing its distance with the middle unstable branch is that the probability function of being about right equilibrium branch is spread. As a result, simulated trace of s and h_e will be more distributed about the up state equilibrium point and sometimes reaches to the area of middle branch and may pass this unstable barrier and switch the whole network from up to down state. This increases s_∞ rapidly but s that is slowly following s_∞, is increased gradually. This is the reason why the simulated trace in Figure 4A reaches to the U-turn tip while passing the unstable barrier, s value in the U-turn tip is a critical value that enough modulatory forces exist in the network; so the network can not be settled in down state any longer and must be switched back to up state. Reaching to this critical value conveys an abrupt increasing of excitatory and inhibitory membrane potential and decreasing of s_∞. After this switching, s is gradually decreased until another jump from up to down state happens.

Open in a separate window

Figure 4

The same description of Figure 3 except that anesthetic drug concentration has been increased to 0.3 mM.

V. Conclusion

In the case where a MFM contains many state variables it is not feasible to study the behavior of the MFM by considering all of its state variables; so rational model order reduction is the first step that must be done. In our EMFM, we adiabatically reduced the model order from nine to three. One of them was a slow variable (s) and others were faster variables (h_e, h_i); so a 3D vector field was necessary for investigation on the role of slow and fast variables. This 3D vector field was reduced to a 2D vector field by combining the two fast variables. We could validate this method by observing that the 2D vector field can be matched with the trajectories of numerically simulated s and h_e signals.

References