| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
INNOVATIVE METHODOLOGY
Department of Biology, Emory University, Atlanta, Georgia
Submitted 27 July 2007; accepted in final form 5 September 2007
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; MacLean et al. 2003; Marder and Goaillard 2006; Swensen and Bean 2005). These maintained functional characteristics could be constancy in the period of a network producing rhythmic output, the time interval between the stimulus onset and the response of a neuron, or many others.
Various factors compensate each other to maintain function. An example of such compensation has been discovered recently in neurons of the lobster stomatogastric ganglion (MacLean et al. 2003; Schulz et al. 2006). In neurons with similar firing properties, there was a linear correlation between the conductances of transient potassium (IA) and hyperpolarization-activated inward (Ih) currents. This correlation accounted for the commonly observed variability of these two currents across preparations. Importantly, the correlation appears to reflect a homeostatic mechanism that regulates IA and Ih. Experimental (MacLean et al. 2003; Swensen and Bean 2005) and modeling (Achard and De Schutter 2006; Goldman et al. 2001; MacLean et al. 2005; Prinz et al. 2004) studies in vertebrates and invertebrates show that multiple combinations of conductances can underlie similar neuronal behavior, suggesting that multiple complex compensatory mechanisms are possible. A traditional approach for finding compensatory mechanisms is based on sensitivity analysis (MacLean et al. 2005; Nygren et al. 1998; Olsen et al. 1995; Paulsen et al. 1982; Weaver et al. 2007). In this analysis, the sensitivity of a specified characteristic to a certain parameter is a ratio of the normalized change in the characteristic to the underlying normalized change in the parameter. Sensitivities can be used to determine approximately how much one parameter must change from its original value to compensate a change of the other parameter. That multiple compensatory covariations exist is well documented (Marder and Goaillard 2006; Rich and Wenner 2007; Ward 2006); what is lacking are systems for determining and understanding these covariations.
In this study, we introduce a general method for describing all putative compensatory covariations of neuronal and network parameters given specified characteristics of neuronal activity to be maintained. The method rigorously specifies exact covariations and linear approximations to them. It is based on a fundamental and powerful mathematical tool, the implicit function theorem. We show that compensatory covariations are possible charts of manifolds composed of parameter sets for which specified characteristics are constant. In particular, we show that to maintain N independent characteristics of neuronal activity affected by a change of one or more parameters, at least N other parameters should change appropriately.
To illustrate the method, we apply it to a model of the smallest functional network of the leech heartbeat central pattern generator (CPG) consisting of two mutually inhibitory neurons (Cymbalyuk et al. 2002; Hill et al. 2001; Olypher et al. 2006; Sorensen et al. 2004) (Fig. 1A). In particular, we show how certain parameters must covary to maintain the period of the network activity and the average spike frequency in bursts. Finally, we discuss applications of our method to the analysis of the homeostatic regulation, model tuning and other problems.
|
, 2007). |
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
of a half-center oscillator from the leech heartbeat CPG, consisting of two identical neurons that mutually inhibit one another. The model, like the living network, produces a periodic pattern of activity where the two neurons burst alternately (Fig. 2). The neurons were modeled according to the Hodgkin-Huxley formalism with five inward and three outward voltage-dependent currents. The inhibitory interactions of the model neurons, as in the living system, have a graded and a spike-mediated component (Angstadt and Calabrese 1991; Ivanov and Calabrese 2000, 2003, 2006a,b; Lu et al. 1997; Olsen and Calabrese 1996). In this study, by the synaptic strength we mean the maximal conductance, 
SynS, of the spike-mediated current. In the examples considered, synaptic strength was always one of the parameters varied. The other parameters were the maximal conductance, h, of a hyperpolarization-activated current, Ih, and a scaling factor, 
, for the inactivation time constant, 
h,CaS, of a slowly inactivating low-threshold calcium current, ICaS
 |
,CaS(V) = 1/(1 + e360(V–0.055)), h,CaS(V) = 0.2 + 5.25/(1 + e–2.50(V+0.043)), V was the membrane potential (in V), t was time (in s) (Hill et al. 2001). By definition, greater (smaller) than one slows down (speeds up) the inactivation. When a parameter of the neuronal model was varied, it was varied in both neurons simultaneously to comply with the simplifying assumption that the two neurons in the half-oscillator model are identical. However, our method can be similarly applied for unilateral variations of parameters.
|
L = 9.9 nS and the reversal potential EL = –63.5 mV instead of L = 8 nS, EL = –60 mV as in the Hill et al. (2001) model (cf. Cymbalyuk et al. 2002). Other parameters of the Hill et al. (2001) model were not changed. The model and the set of chosen parameters are called canonical here.
Another change in the Hill et al. (2001) model was aimed at making it smooth with respect to the variables and parameters. By smoothness here and in the following text, we mean at least first-order differentiability. In other words, the right-hand sides of the model system must change continuously with changes of the variables and parameters, and the velocity of these changes must also be continuous. We consider this requirement as quite natural and meeting physiological intuition. The original model has only one place that required a change to meet this requirement—the expression for the total calcium current ICa = max(0, –ICaF – ICaS – A), where ICaF and ICaS are the rapidly inactivating and slowly inactivating slow low-threshold calcium currents, respectively, and A is a dynamic threshold (see details in Hill et al. 2001). To make the system smooth, this original nonsmooth function for ICa was replaced here by a smooth sigmoid function 
This substitution, as expected, had negligible effect on the solutions of the system.
To quantify the system's activity for a particular set of parameters, the model was simulated for 110 s. Measured characteristics were burst period, T, defined as in Hill et al. (2001) as the time between the middle spike of one burst and the middle spike of the next burst, and average intraburst spike frequency, F, defined as a number of spikes in a burst divided by the period. The period was calculated with the middle spikes rather than the first or last spikes because it yielded the smallest coefficient of variation. These two characteristics were chosen because in the living system, the period of these oscillations set heartbeat period, and the intraburst spike frequency determines the level of inhibition both within the half-center oscillators and between the oscillator interneurons and their motor neuron targets (Hill et al. 2001). The first 30 s of the simulation was considered as an interval sufficient for stabilizing the pattern and were discarded from analysis. Convergence to a stable regime was confirmed by low variation of the period over the next 80 s of simulation. The initial conditions for all simulations were the same and taken from the stable regime corresponding to the canonical parameters. Simulations were performed with the Matlab (MathWorks, Natick, MA) solver ode15s with the absolute and relative tolerances 10–9 and 10–8, respectively. Matlab code is available from http://calabreselx.biology.emory.edu/pub/HC.zip.
Partial derivatives of activity characteristics with respect to parameters were first assessed with standard numeric formulas: by varying a parameter and finding the ratio of the characteristic's change over the parameter's change, repeating this for doubled variations of the parameters, and using the Richardson extrapolation to improve the accuracy of the result (Press et al. 1993). We also estimated the partial derivatives by calculating the characteristics for the original values of the parameters and ±8, ±6, ±4, and ±2% deviations from the original values and then applying a linear fit to these values of the characteristics. Both methods gave similar estimates of the partial derivatives. Complication arose when considering the partial derivative of the average intraburst spike frequency, F, with respect to SynS. The changes of F for 
8% changes of CaS were too small to obtain a reliable linear fit. Thus we used a linear fit of the values of F corresponding to ±20 and 0% changes of SynS to estimate the derivative and obtained a number with the absolute value less then 0.005 (see following text).
Implicit function theorem
Our main analytical tool was the implicit function theorem (Rudin 1976). The theorem states that under certain conditions, equalities involving several variables imply functional dependencies between the variables (APPENDIX A). For example, if the conditions of the theorem are met, an equality C(p1, p2) = 0, true for particular values p1 = p1*, p2 = p2*, implies that near the point (p1*, p2*) there exists a unique function p2= f(p1) such that p2* = f(p1*), C(p1, f(p1)) = 0, and 
The first condition of the theorem is the differentiability of characteristics with respect to parameters, which in the example considered means that partial derivatives 
C(p1*, p2*)/p1, and C(p1*, p2*)/p2 must exist. The second condition of the theorem is the nondegenerate dependency of characteristics with respect to parameters that are supposed to be dependent on other parameters. In the current example, this condition is fulfilled if C(p1*, p2*)/p2 
0.
Algorithm for specifying the variations of model parameters that maintain activity characteristics
The following algorithm is a general one. It can be applied to an arbitrary neuronal or network model to determine coordinated changes of n parameters that do not affect m characteristics. The main assumption of the algorithm is a smooth dependence of characteristics on parameters.
First, select parameters p1, p2, ..., pn of the model to be varied near some point p* = (p1*, p2*, ..., pn*) in the parameter space.
Second, select a set of activity characteristics C(p) = (C1(p), C2(p), ..., Cm(p)) to be maintained. By subtracting target values from Ci, we can assume that C(p*) = 0, and the goal is to find p near p* such that C(p) = 0. Let partial derivatives Ci(p*)/pj, i = 1, ..., m, j = 1, ..., n exist.
Third, assume that for y = (p1, ..., pm), the determinant
 |
y = –[Cy'(p*)]–1·Cx'(p*)·x, where y = (p1 – p1*, ..., pm – pm*), x = (pm+1 – pm+1*, ..., pn –pn*), and
 |
Definitions: we shall call p1, ..., pm compensating parameters, pm+1, ..., pn compensated parameters, and pi = fi(pm+1, ..., pn), 1 i m, compensatory functions.
Notes on the algorithm
The calculation of the gradients of the selected characteristics is a key computational step. This step is in common with multiple computational procedures, such as for example the steepest descent method. The difference here is that we are not interested in the direction determined by the gradient, but in the direction, or directions in a multidimensional space, orthogonal to the gradient.
The condition in step 3 of the algorithm means that there is no characteristic for which the gradient in p* is a linear combination of the gradients of other characteristics. This condition does not hold if there is a linear dependence among characteristics, in the simplest case when Cj(p) = k·Cl(p) for some 1 j, l m, and k 0. Only independent characteristics can be used for the analysis. Another case where the condition does not hold is when there is a pair of linearly correlated parameters, e.g., p1 = 
·p2 for some 0. Only one parameter from such a pair must be kept in the set of compensating parameters.
The condition in step 3 may not be satisfied at critical (bifurcational) points of the parameter space. Indeed at such points, a solution to the system, for which the dynamics (characteristics) are specified, either disappears or changes abruptly. These points limit areas where functional dependencies between parameters that maintain specified characteristics can be defined with the help of the implicit function theorem. For the period of the activity, a corresponding result is well known (APPENDIX B).
Step 5 does not specify a method for refining the isomanifold. In this study, we used standard Matlab routines. To find the isomanifold, on which the period T of the model is equal to the target value T*, we used the function fzero. fzero finds a zero of a function, in our case T – T*, by a combination of bisection, secant, and inverse quadratic interpolation methods. To find the isomanifold, on which the period T of the model is equal to the target value T*, and the spike frequency F is equal to the target value F*, we used the function fminsearch. fminsearch minimizes a function, in our case (T – T*)2 + (F – F*)2, by the simplex search method.
The implicit function theorem provides a clear geometric description of the set of parameter combinations for which specified characteristics of activity are constant. The theorem shows that these sets form smooth, n-m dimensional manifolds in the parameter space. These manifolds are intersections of m (n – 1)-dimensional manifolds each of which is a level set of a particular characteristic. Because the latter manifolds can be called isosurfaces, we'll refer to their intersection as an isomanifold. All possible compensatory covariations of parameters, with some m parameters being functions of the other n-m parameters, are possible charts (Arnold 1978) of this isomanifold. These charts are well defined according to the implicit function theorem.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
To illustrate our method for determining compensations of functional characteristics through parameter variation, we consider examples from our own work. In a leech half-center oscillator from the heartbeat CPG, the strength of the mutually inhibitory synapses varies over a more than twofold range yet cycle period, spike frequency, and other functionally important output characteristics can remain very similar (Tobin and Calabrese 2005) (Fig. 1). Similar variations of the synaptic strength caused significant changes (10%) in the period of our canonical model of this half-center oscillator (Fig. 2). We later consider another functionally important characteristic—the average spike frequency in a burst. Spike frequency determines the accumulated impact of the spike-mediated inhibition in the half-center oscillator (Hill et al. 2001). In the model, corresponding changes in synaptic strength also changed the spike frequency but less than the period and nonmonotonically (Fig. 2A, right).
What parameters of the model can co-vary with SynS in such a way that the period will remain constant, as in the experiments of Tobin and Calabrese (2005)? Multiple parameters, in addition to SynS, strongly affect the period in the model (APPENDIX C, Table 1). Among these parameters are the maximal conductances of the low-threshold slowly inactivating calcium current, persistent potassium current and other currents (Cymbalyuk et al. 2002; Hill et al. 2001; Olypher et al. 2006; Sorensen et al. 2004). According to the general picture described in METHODS, variations of all these parameters compensating variations of SynS are possible parameterizations of the manifold in the parameter space where the period is constant. For the current analysis, we chose the maximal conductance, h, of the hyperpolarization-activated current; any other parameter can be analyzed similarly. Next, we just followed the steps of the general algorithm.
|
h, p2 = SynS; p1* = 4 nS, p2* = 150 nS.
Second, a specified characteristic was C = T(p) – T(p*), where T was the cycle period, T(p*) = 7.91 s, C(p*) = 0. A smooth dependence of T, and hence C, from SynS and h follows from the theorem, described in METHODS and illustrated by Fig. 3A.
|
h, x = SynS. Then det[Cy'(p*)] = T(p*)/h = –0.63 s/nS 0 (Fig. 3A).
Fourth, according to the implicit function theorem, near SynS* = 150 nS and h* = 4 nS there exists a unique function h = f(SynS) such that for each pair (f(SynS), SynS), the cycle period T is the same as at the point SynS = SynS*, h = h*. The function f is smooth and has a linear approximation (the tangent)
 | (1) |
= –[T(p*)/h]–1·[T(p*)/SynS] = –0.077/(–0.63) = 0.012.
Fifth, we used linear approximation (Eq. 1) to calculate the exact compensatory functional dependency f between h and SynS (Fig. 3B) near SynS* = 150 nS, and h* = 4 nS. For larger values of SynS, we used the parabolic approximation of the function h = f(SynS) based on its values near SynS* = 150 nS, and h* = 4 nS. Figure 4A shows that the linear approximation (Eq. 1) provided compensation of the period within 2% over a wide range of varied SynS, which is a consequence of the fair approximation of f by its tangent over the same range. When h was determined by the functional dependency shown in Fig. 3B, the compensation of the period was nearly perfect (Fig. 4B1). Both approximate and exact compensatory covariations of h and SynS affected the spike frequency much more than corresponding variation of SynS alone (Fig. 4, A2 and B2).
|
SynS and h have their canonical values. Near the point (h*, SynS*), this section is a smooth one-dimensional manifold, which is simply a smooth curve (cf. Fig. 3B). Obviously, we could parameterize this curve not by SynS but by h, and interpret the resulting function, inverse to f, as one determining variations of SynS compensating variations of h. Maintaining both period and spike frequency of the network in the face of parameter variation
According to the theory developed, to maintain both the period, T, and the spike frequency, F, when SynS is varied, compensatory covariations of at least two parameters are necessary and sufficient. As one of the many potential pairs of parameters, we chose h and a scaling factor, , of the time constant of the slow low-threshold calcium current inactivation, hCaS (see METHODS). As we have shown recently, strongly correlates with the burst duration of the neuron (Olypher et al. 2006). In this case, the application of the general algorithm was as follows.
First, selected parameters were p = (p1, p2, p3) with p1 = h, p2 = , p3 = SynS; p* = (p1*, p2*, p3*), p1* = 4 nS, p2* = 1, p3* = 150 nS. Note that p2 is dimensionless, whereas p1, p3 are in nanoSiemens. It is easy to check that the algorithm is not sensitive to differences in units among parameters and/or characteristics.
Second, specified characteristics were C1 = T(p) – T(p*), where T(p*) = 7.91 s, and C2 = F(p) – F(p*), where F(p*) = 8.82 Hz. As earlier, a smooth dependence of T, and hence C, from SynS, h, and follows from the general theory. With regards to the frequency, we approximated it by a smooth function (Fig. 5A2).
|
h, y2 = , x = SynS. Then
 |
h = 
(SynS), = 
(SynS) such that for each triple [(SynS), (SynS), SynS], the cycle period T, and the spike frequency F are the same as at p*. The functions and are smooth. Their linear approximations follow from the system
 |
 | (2) |
SynS*, 1.5 SynS*]. There were no difficulties in refining this approximation for the values of SynS close to 150 nS (Fig. 6). To calculate the points on the isomanifold for the values of SynS >250 nS, we used linear approximations based on points already found on the manifold. Note that the compensating dependency h = h(SynS) is not monotonic (Fig. 6B) as in the case when only the period was maintained (Fig. 3B). This nonmonotonicity, to be considered in detail elsewhere, is in contrast to the monotonic and almost linear compensating dependency = (SynS) (Fig. 6C). Figure 7 shows that the activity in the compensated model is very close to the activity of the canonical model.
|
|
SynS, h, and had their canonical values, was a smooth curve. We parameterized this curve by SynS and interpreted the resulting parameterization as variations of h and compensating variations of SynS. Equation 2 is a linear approximation to that parameterization. We could of course consider h as a free parameter and parameterize the same smooth curve with h. That would give us variations of SynS and compensating variations of h.
To explore the individual contributions of h and to the approximate, linear compensation, we first made simulations with h varied according to h = 0.00049·SynS and fixed h. Figure 5B shows that h compensated both the period and the spike frequency. Indeed, from the partial derivatives T(p*)/h = –0.63, and F(p*)/h, it followed, that 1% increase of h should lead to 0.32% decrease and 0.39% increase of T and F, respectively. Next we made simulations with fixed h = 4 nS and varying according to = –0.00080·SynS (Fig. 5C). Corresponding partial derivatives implied a notably greater effect of on T (2.9% increase) than on F (0.49% increase) of F with 1% increase of . For small deviations of SynS from its canonical value, the contributions of h and are approximately independent and proportional to changes in SynS, showing that the compensation works in linear approximation.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Our approach requires a smooth dependency of activity characteristics on parameters. Such smoothness is often intrinsic to neuronal systems when they are in stable functional states. The conclusions about how parameters compensate each other, developed in this study, can thus be used even without regard to the specific mathematical model describing a particular neuron or neuronal network.
The application of the implicit function theorem to finding level sets of a characteristic of neuronal activity is not new. In particular, it is the basis of the parameter continuation technique implemented in dynamical systems software packages, like XPP-AUTO (Ermentrout 2002), and CONTENT (Kuznetsov et al. 1996). The theorem, when its conditions are met, guarantees the existence and uniqueness of a continuation and provides a linear approximation to compensatory functions which can be refined. In practice, parameter continuation is mostly used for finding isocurves, i.e., level sets of a characteristic in a plane of two parameters. Here, we consider a general case with multiple characteristics and parameters, focus on the geometrical picture provided by the theorem, and most importantly, connect the theory with the general problem of neuronal variability.
Predictions
Our general approach allows us to formulate several predictions related to the analysis of the homeostatic regulation, model tuning, and data analysis. Specifically, if a particular homeostatic mechanism maintains m independent characteristics of neuronal activity, then our approach predicts that at least m parameters must be changed as a response to perturbation in one or more other parameters of the system. Correspondingly, if during model tuning, one parameter is updated, then to return to target values of m characteristics of the model activity, at least m other parameters have to be readjusted.
Our study shows that correlations between pairs of parameters in a database of simulated models with similar behavior are generically rare but nevertheless possible. They are rare because a correlation between two parameters adds an extra constraint on the set of parameter combinations for which characteristics remain constant. Geometrically, this constraint means that an isomanifold on which characteristics are constant must belong to a hypersurface orthogonal to the hyperplane of the two parameters. On the other hand, we have demonstrated an almost perfect linear correlation between the two parameters and SynS in the half-center oscillator model explored here (Fig. 6C) under the condition that the burst period and spike frequency in the model remain constant with variations of these two parameters and h.
Our theory also implies that a biological control mechanism that coregulates (linear correlation) different ion channels (conductance parameters), regardless of their functional relation to one another, may effectively limit the number of activity characteristics that can be maintained. Effectively the parameters become nonindependent in regulating activity. In neurons of the lobster stomatogastric ganglion, linear correlation between the conductances of IA and Ih have been observed (MacLean et al. 2003; Schulz et al. 2006). Overexpression of IA by mRNA injection leads to corresponding increases in Ih expression and relative constancy of cell activity (burst period and spike frequency). One is tempted to conclude from our analysis that such correlation means that these currents exclusively control burst activity. We would argue that rather than indicating that other parameters do not influence burst activity, it is more likely that the epigenetic mechanism was selected by evolution as more efficient than other potential compensatory mechanisms involving other parameters. At the same time such a mechanism can limit the system's flexibility. Perhaps such epigenetic homeostatic mechanisms are useful, however, for setting baseline activity. Modulation may break such linkages by affecting one current more than the other either by channel modification or by directly affecting gene expression, and thus alter this homeostatically set baseline activity (see, e.g., Khorkova and Golowasch 2007). Moreover, combinations of such epigenetic mechanisms that effect different correlations between the two parameters may restore flexibility (see, e.g., Brezina et al. 1996).
Relation to previous work
Our approach extends sensitivity analysis (MacLean et al. 2005; Nygren et al. 1998; Olsen et al. 1995; Paulsen et al. 1982; Weaver et al. 2007), compliments the "brute force" approach (Achard and De Schutter 2006; Goldman et al. 2001; MacLean et al. 2005; Prinz et al. 2004), and is "orthogonal" to bifurcational analysis (Izhikevich 2006). In sensitivity analysis, sensitivities determine approximately how much one parameter should change from its original value to compensate a deviation of another parameter from its original value. Sensitivities approximate partial derivatives of a characteristic at a point of interest in the parameter space. When this approximation is accurate, sensitivity analysis provides a linear approximation to actual compensatory covariations of parameters. Traditionally, the maintenance of only one activity characteristic by compensatory covariations of a pair of parameters has been considered. Here we consider multiple characteristics and parameters and show how to reveal exact compensatory functional dependencies between parameters.
The brute force approach was critical in revealing the fact that neuronal and network models with multiple sets of parameter can have similar behavior (Goldman et al. 2001; Prinz et al. 2004). In particular, the sampling of the large volume of parameter space proved that these parameter sets can be very different. Our approach is local. We start with one point in the parameter space and explore its vicinity. The payoff is that we get complete information on these parameter sets near the point. Our simulations, described here, show that the vicinity in which our analysis is applicable can be reasonably large (cf. Fig. 4B1). Furthermore, the isomanifolds the existence and uniqueness of which are established locally can be extended out of that vicinity.
Sampling the parameter space of models has also been used for analyses of the structure of parameter sets for which a system's behavior is the same. For example, Goldman et al. (2001), MacLean et al. (2005), and Taylor et al. (2006) found that these sets have often a layered structure. This observation is in accord with the general picture developed here. Indeed, manifolds corresponding to close characteristic values should not intersect because the intersection would mean that for the same set of parameters, activity characteristics would be different. From the layered structure, Goldman et al. (2001) deduced that there are sensitive and insensitive directions in the parameter space. In the sensitive direction, a characteristic changes most of all, whereas in the insensitive direction, it changes least of all. Our analysis puts these observations into a rigorous mathematical framework. Namely, if the smooth manifold of parameters is one dimensional, then at every point it has a gradient, a vector with components that are partial derivatives of the characteristic with respect to parameters. By definition, the gradient shows a direction, in which the characteristic changes most of all, or the most sensitive direction. Directions, orthogonal to the gradient, are obviously insensitive directions.
Achard and De Schutter (2006), who simulated a multicompartmental model, showed that once there is a point in the parameter space with values close to the target values of the specified characteristics then in the vicinity of this point there are other points where the model has characteristics that are also close to those target values. From our analysis, it follows that near a generic point there must be a whole manifold of parameters for which the characteristics are the same.
Our approach is "orthogonal" to bifurcational analysis, which is about small parameter changes leading to significant changes in dynamics (activity). The focus of our approach is to seek variations of parameters that have no or small effect on model dynamics. Given that characteristics of the dynamics change drastically when a system experiences a transition to another state, linear approximations to compensatory covariations of parameters work better the further parameters are from such critical values. Figure 8A illustrates such a transition from the alternating bursting to spiking in the model considered here.
|
To illustrate the method, we applied it to a model of the smallest functional network of the leech heartbeat CPG consisting of two mutually inhibitory neurons. Motivated by the experimental finding by Tobin and Calabrese (2005) (Fig. 1), we varied the maximal conductance, SynS, of the spike-mediated inhibition between the neurons and explored how other parameters might compensate variations of SynS. Our choice of the maximal conductance of the hyperpolarization-activated current, h, and the scaling factor, , of the inactivation time constant of the slow low-threshold calcium current as potential compensating parameters was influenced by our recent studies of the role of these parameters in shaping pattern of the activity in this particular network (Olypher et al. 2006; Sorensen et al. 2004) and the analysis of covariations of h and SynS in hybrid half-center oscillators formed from stomatogastric neurons (Grashow and Marder 2006).
Although extensive model simulations and physiological intuition strongly suggested that these three parameters might change in a coordinated way to maintain the burst period and spike frequency, it was not clear that such maintenance was possible in the model, given the complex nonlinear interactions between the parameters and characteristics. The theory presented in this study firmly establishes the existence and uniqueness of such compensatory functional dependencies.
We generated here a testable prediction on concordant changes of h, SynS, and . This prediction can be verified in a hybrid system consisting of a living, pharmacologically isolated, heart interneuron connected with artificial synapses to a model heart interneuron running in real-time (Olypher et al. 2006; Sorensen et al. 2004). In this setting, the endogenous hyperpolatization-activated current and calcium current of the living interneuron can be pharmacologically blocked and substituted by corresponding artificial conductances with regulated h, SynS, and .
Applicability and perspectives of the method
In the method developed here, characteristics and the parameters maintaining them can have arbitrary nature with regards to units and origin. In particular, a close analysis of the algorithm shows that at every step only terms with similar units are summed. Both parameters and characteristics can originate in ordinary differential equations, partial differential and integral equations, mappings, multicompartmental models, or continuous medium network models. The main requirement is that the characteristics considered must smoothly depend on parameters. However, even this requirement can be lifted, if not the characteristic itself, but its smooth approximation is considered and if this approximation depends on parameters smoothly. For example, in the interval [0.5SynS, 1.5SynS], the number of spikes in the noncompensated model explored here increased from 33 to 36 meaning that there were two points where the partial derivative of the spike frequency with respect to SynS lost its continuity. However, in that interval the compensated spike frequency was <2% different from the canonical value indicating that the compensation worked.
The accuracy of the linear approximation of the isomanifold depends on the accuracy of the estimates of partial derivatives of characteristics with respect to parameters. Figure 8B illustrates a case where the period in the model for a particular choice of the leak current parameters is so variable that a reliable estimate of the period's partial derivative with respect to the reversal potential of the leak current, EL, is impossible.
Accurate estimates of partial derivatives of characteristics do not guarantee that specific parameters can be considered as compensated or compensating. The critical condition is formulated in the step 3 of the algorithm. If the corresponding determinant is equal to zero within the numeric accuracy employed, then the set of compensated and compensating parameters must be rearranged/reduced.
A particular numeric method for finding an isomanifold must be chosen depending on the problem at hand. For example, for finding isosurfaces there are effective algorithms developed for data visualization, such as the Marching Cubes algorithm and its generalizations (Bhaniramka et al. 2004). In the case of multiple characteristics, the intersection of the isosurfaces must be considered. For this case, there may be specific and effective algorithms as well; however, we are not aware of such algorithms.
In this study, we used the implicit function theorem for revealing functional relations between finite-dimensional vectors of parameters. The theorem has a more general form, defined for functional spaces (Wouk 1979) with the premises and statements of the theorem reformulated correspondingly. The theorem states the local existence and uniqueness of an operator (nonlinear) for compensatory changes not of parameters but of functions affecting specified (functional) characteristics. In theoretical analyses, this general form of the implicit function can be useful in proving existence and uniqueness of implicitly defined connections between various functions affecting neuronal and network activity.
|
APPENDIX A: IMPLICIT FUNCTION THEOREM |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Theorem (Rudin 1976). Consider a mapping F:U
Rn defined in a neighborhood U of point (x0, y0) 
Rn+m, where m and n are natural numbers. Let F be continuously differentiable, F C(p)(U), F(x0, y0) = 0, the Jacobian of F with respect to y is an invertible matrix, det[Fy'(x0, y0)] 0. Then there exists the m + n dimensional box I = Ixm x Iyn 
U, where Ixm = {x Rm:|x – x0| < }, Iyn = {y Rn:|y – y0| < β}, and and β are some positive numbers, and such a mapping F C(p)(Ixm x Iyn) that F(x, y) = 0 
y = f(x) for each point (x, y) Ixm x Iyn, and f'(x0) = –[Fy'(x0, y0)]–1[Fx'(x0, y0)].
Here, C(p)(U) is a standard notation for the space of functions, which are continuous and have continuous partial derivatives up to order p in U; Fy'(x0, y0), and Fx'(x0, y0) are matrices made of partial derivatives of the vector function F with respect to components of vectors y and x in (x0, y0).
|
APPENDIX B: SMOOTHNESS OF THE PERIOD |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
). Let x, F be real vectors. Let F (x,µ) be continuous for small |µ| and for x on some d –dimensional domain and have continuous partial derivatives with respect to the components of x. When µ = 0, let x' = F (x, µ) have a nonconstant solution x = g0(t) of period p0 > 0 such that the linearized system y' = A(t)y, where A(t) 
xF [g0(t), 0], has exactly one characteristic root equal to 1. Then, for small |µ|, the system x' = F (x, µ) has a unique periodic solution x = x (t, µ) with a period p(µ), depending on µ, such that x(t, µ) is near g0 (t) and the period p(µ) is near p0; furthermore x(t, µ), p(µ) are continuous, x(t, 0) = g0(t), and p(0) = p0.
Further analysis shows that when F (x, µ) is differentiable, x(t, µ), p(µ) are also differentiable (Hartman 1964).
|
APPENDIX C: PERIOD SENSITIVITIES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
GRANTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: A. V. Olypher, Dept. of Biology, Emory University, 1510 Clifton Rd. N.E., Atlanta, GA 30322 (E-mail: aolypher{at}biology.emory.edu)
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIX A: IMPLICIT FUNCTION...APPENDIX B: SMOOTHNESS OF...APPENDIX C: PERIOD SENSITIVITIESGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Angstadt JD, Calabrese RL. Calcium currents and graded synaptic transmission between heart interneurons of the leech. J Neurosci 11: 746–759, 1991.[Abstract]
Arnold VI. Ordinary Differential Equations. Cambridge, MA: MIT Press, 1978.
Bhaniramka P, Wenger R, Crawfis R. Isosurface construction in any dimension using convex hulls. IEEE T Vis Comput Gr 10: 130–141, 2004.[CrossRef]
Brezina V, Orekhova IV, Weiss KR. Functional uncoupling of linked neurotransmitter effects by combinatorial convergence. Science 273: 806–810, 1996.[Abstract]
Buzsaki G, Csicsvari J, Dragoi G, Harris K, Henze D, Hirase H. Homeostatic maintenance of neuronal excitability by burst discharges in vivo. Cereb Cortex 12: 893–899, 2002.
Cymbalyuk GS, Gaudry Q, Masino MA, Calabrese RL. Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J Neurosci 22: 10580–10592, 2002.
Ermentrout E. Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students. Philadelphia: SIAM, 2002.
Goldman MS, Golowasch J, Marder E, Abbott LF. Global structure, robustness, and modulation of neuronal models. J Neurosci 21: 5229–5238, 2001.
Grashow R, Marder E. Assaying intrinsic property variability in biological and model stomatogastric neurons using the dynamic clamp. Soc. Neurosci Abstr 649.18/W27, 2006.
Hartman P. Ordinary Differential Equations. New York: Wiley, 1964.
Hill AA, Lu J, Masino MA, Olsen OH, Calabrese RL. A model of a segmental oscillator in the leech heartbeat neuronal network. J Comput Neurosci 10: 281–302, 2001.[CrossRef][Web of Science][Medline]
Ivanov AI, Calabrese RL. Intracellular Ca2+ dynamics during spontaneous and evoked activity of leech heart interneurons: low-threshold Ca currents and graded synaptic transmission. J Neurosci 20: 4930–4943, 2000.
Ivanov AI, Calabrese RL. Modulation of spike-mediated synaptic transmission by presynaptic background Ca2+ in leech heart interneurons. J Neurosci 23: 1206–1218, 2003.
Ivanov AI, Calabrese RL. Graded inhibitory synaptic transmission between leech interneurons: assessing the roles of two kinetically distinct low-threshold Ca currents. J Neurophysiol 96: 218–234, 2006a.
Ivanov AI, Calabrese RL. Spike-mediated and graded inhibitory synaptic transmission between leech interneurons: evidence for shared release sites. J Neurophysiol 96: 235–251, 2006b.
Izhikevich EM. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting (Computational Neuroscience). Cambridge, MA: MIT Press, 2006.
Khorkova O, Golowasch J. Neuromodulators, not activity, control coordinated expression of ionic currents. J Neurosci 27: 8709–8718, 2007.
Kuznetsov YA, Levitin VV, Skovoroda AR. Continuation of Stationary Solutions to Evolution Problems in CONTENT. Report AMR9611. Amsterdam: Centrum/Voor Wiskunde en Informatica, 1996.
Lu J, Dalton JFt, Stokes DR, Calabrese RL. Functional role of Ca2+ currents in graded and spike-mediated synaptic transmission between leech heart interneurons. J Neurophysiol 77: 1779–1794, 1997.
MacLean JN, Zhang Y, Goeritz ML, Casey R, Oliva R, Guckenheimer J, Harris-Warrick RM. Activity-independent coregulation of IA and Ih in rhythmically active neurons. J Neurophysiol 94: 3601–3617, 2005.
MacLean JN, Zhang Y, Johnson BR, Harris-Warrick RM. Activity-independent homeostasis in rhythmically active neurons. Neuron 37: 109–120, 2003.[CrossRef][Web of Science][Medline]
Marder E, Goaillard JM. Variability, compensation and homeostasis in neuron and network function. Nat Rev Neurosci 7: 563–574, 2006.[CrossRef][Web of Science][Medline]
Nygren A, Fiset C, Firek L, Clark JW, Lindblad DS, Clark RB, Giles WR. Mathematical model of an adult human atrial cell: the role of K+ currents in repolarization. Circ Res 82: 63–81, 1998.
Olsen OH, Calabrese RL. Activation of intrinsic and synaptic currents in leech heart interneurons by realistic waveforms. J Neurosci 16: 4958–4970, 1996.
Olsen OH, Nadim F, Calabrese RL. Modeling the leech heartbeat elemental oscillator. II. Exploring the parameter space. J Comput Neurosci 2: 237–257, 1995.[CrossRef][Web of Science][Medline]
Olypher AV, Calabrese RL. Local variations of parameters in the model half-center oscillator maintaining constant period. Soc Neurosci Abstr 649.16/W25, 2006.
Olypher AV, Calabrese RL. Variations of neuronal parameters that do not change network output. Proceedings of 16th Annual Computational Neuroscience Meeting, Toronto, Ontario, p. 237, 2007.
Olypher AV, Cymbalyuk G, Calabrese RL. Hybrid systems analysis of the control of burst duration by low-voltage-activated calcium current in leech heart interneurons. J Neurophysiol 96: 2857–2867, 2006.
Paulsen RA Jr, Clark JW Jr, Murphy PH, Burdine JA. Sensitivity analysis and improved identification of a systemic arterial model. IEEE Trans Biomed Eng 29: 164–177, 1982.[Web of Science][Medline]
Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical Recipes in C. The Art of Scientific Computing (ed. 2). Cambridge, UK: Cambridge Univ. Press, 1993.
Prinz AA, Bucher D, Marder E. Similar network activity from disparate circuit parameters. Nat Neurosci 7: 1345–1352, 2004.[CrossRef][Web of Science][Medline]
Rich MM, Wenner P. Sensing and expressing homeostatic synaptic plasticity. Trends Neurosci 30: 119–125, 2007.[CrossRef][Web of Science][Medline]
Rudin W. Principles of Mathematical Analysis. New York: McGraw-Hill, 1976.
Schulz DJ, Goaillard JM, Marder E. Variable channel expression in identified single and electrically coupled neurons in different animals. Nat Neurosci 9: 356–362, 2006.[CrossRef][Web of Science][Medline]
Sorensen M, DeWeerth S, Cymbalyuk G, Calabrese RL. Using a hybrid neural system to reveal regulation of neuronal network activity by an intrinsic current. J Neurosci 24: 5427–5438, 2004.
Swensen AM, Bean BP. Robustness of burst firing in dissociated Purkinje neurons with acute or long-term reductions in sodium conductance. J Neurophysiol 25: 3509–3520, 2005.
Taylor AL, Hickey TJ, Prinz AA, Marder E. Structure and visualization of high-dimensional conductance spaces. J Neurophysiol 96: 891–905, 2006.
Tobin A-E, Calabrese RL. Myomodulin increases Ih and inhibits the Na/K pump to modulate bursting in leech heart interneurons. J Neurophysiol 94: 3938–3950, 2005.
Ward NS. Compensatory mechanisms in the aging motor system. Ageing Res Rev 5: 239–254, 2006.[CrossRef][Web of Science][Medline]
Weaver CM, Gamkrelidze G, Baker R, Wearne SL. Sensitivity analysis enables comparison of how realistic morphology and other intrinsic properties influence neuronal firing. Proceedings of 16th Annual Computational Neuroscience Meeting, Toronto, Ontario, p. 65, 2007.
Wouk A. A Course of Applied Functional Analysis. New York: Wiley, 1979.
This article has been cited by other articles:
|
|
 |
|
  A. L. Taylor, J.-M. Goaillard, and E. Marder How Multiple Conductances Determine Electrophysiological Properties in a Multicompartment Model J. Neurosci., April 29, 2009; 29(17): 5573 - 5586. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
C. Gunay, J. R. Edgerton, and D. Jaeger Channel Density Distributions Explain Spiking Variability in the Globus Pallidus: A Combined Physiology and Computer Simulation Database Approach J. Neurosci., July 23, 2008; 28(30): 7476 - 7491. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
