| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Department of Neurobiology and Anatomy, W.M. Keck Center for the Neurobiology of Learning and Memory, The University of Texas Medical School at Houston, Houston, Texas
Correspondence: Address reprint requests to John H. Byrne, Tel.: 713-500-5602; E-mail: john.h.byrne{at}uth.tmc.edu.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
–3). One of the key goals of systems biology is to develop tools that may lead to a better understanding of nonlinear behaviors of biochemical and gene networks at both the molecular and systems levels (4).
A gene and protein network can be described by a mathematical model consisting of ordinary differential equations (ODEs). Bifurcation analysis is a mathematical technique that enables determination of the stability of a system with respect to a parameter (5,6). Bifurcation diagrams describe the dependence of a state variable on a continuous change in a chosen system parameter, termed a bifurcation parameter. A bifurcation is said to take place when there is a change in the number or the stability of solutions of a system. For example, steady-state solutions for the values of the dependent variables may appear, disappear, change stability, or multiple steady-state solutions may coexist. The coexistence of multiple steady-state solutions at a particular value of a bifurcation parameter is termed multiplicity (or multistability, if the solutions are stable to small perturbations). Multiplicity may occur with supralinear and positive feedback interactions among components of a system of coupled ODEs, and oscillatory dynamics may be sustained if a negative feedback loop of interactions is present. Singularity theory (for review see (7,8)) provides a systematic framework to determine how many topologically distinct bifurcation diagrams exist in a nonlinear dynamic system, and to partition the multidimensional parameter space of the model into regions in which different types of bifurcation diagrams exist. This information can be used to classify control parameters, which play a crucial role in determining system dynamics by governing transitions between qualitatively different bifurcation diagrams.
This study applies bifurcation and singularity analysis to a relatively complex signal transduction and gene network that underlies the induction of long-term memory (LTM) to examine model dynamics and determine control parameters. Sensorimotor neuron synapses of the mollusk Aplysia have been utilized extensively as a model system for the study of the cellular and molecular processes underlying learning and memory (9–13). These synapses exhibit both short- and long-term facilitation after exposure to 5-HT. Long-term facilitation (LTF) requires both activation of protein kinase A (PKA) and transcription. Molecular processes that underlie LTF have been studied in detail (for review see (13)). LTF is a correlate of long-term sensitization (LTS) of Aplysia defensive withdrawal reflexes, a form of long-term memory (LTM) (13–16). Pettigrew et al. (17) developed a mathematical model of biochemical processes that underlie the induction of LTF. The model simulates experimental observations of the time courses of activity of protein kinases essential for LTF, specifically PKA and the MAP kinase isoform termed extracellular-regulated kinase (ERK), after temporally spaced or massed exposures to 5-HT. The model also simulates induction of a gene product (Aplysia ubiquitin hydrolase, or Ap-Uch) that is essential for LTF.
This study addresses the following questions regarding the nonlinear features of the Pettigrew model:
In analyzing the dynamics of the model, this study focuses on the role of a positive feedback loop involving reciprocal induction of PKA activity and Ap-uch gene expression. The results indicate that the positive feedback loop can generate a robust and irreversible switch between a basal steady state of low PKA activity and gene expression, and a stimulated steady state of high PKA activity and gene expression. A prolonged steady state of high gene expression might contribute to the strengthening of synapses between neurons, and therefore to the induction and maintenance of LTM. Singularity theory was used to construct classification diagrams of the parameter space. These diagrams help identify molecular processes that might be manipulated, perhaps pharmacologically, to enhance the induction and maintenance of LTM. These methods of systematic analysis will also be helpful to understand a broad class of complex biological phenomena amenable to modeling, such as development and cellular differentiation (18–23), or the cell cycle (3).
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
) was used (see Table 1). The Pettigrew model (17) was developed to study biophysical mechanisms that may be implicated in learning and memory. Analysis of this model shows that it does not exhibit multiplicity of solutions (the coexistence of multiple steady-state solutions) with the parameter values previously published in Pettigrew et al. (17). A goal of this study is to illustrate a methodology for finding various types of bifurcation diagrams and their location in the parameter space. By altering four parameters, i.e., kApSyn, kApSynBasal, Kpka, and Kerk (new values given in Table 1), the model exhibits multiplicity and, therefore, a more diverse set of bifurcation diagrams. With these four new parameter values, the model continues to respond to stimuli as in Pettigrew et al. (17).
|
).) Pettigrew et al. (17) assumed specific molecular identities for the two transcription factors. TF-1 was assumed to be the cAMP response element binding protein (CREB1), and TF-2 was assumed to be CREB2. However, recent results (24) suggest that induction of Ap-uch by 5-HT may not be mediated directly by CREB1 and CREB2, but rather by as yet unidentified transcription factors. Nevertheless, Ap-uch induction is regulated via the cAMP signaling pathway (25), and [cAMP] elevation activates both PKA and ERK in neurons (26,27). Therefore, it remains plausible that PKA and ERK induce Ap-uch by phosphorylating transcription factors, as the model assumes (Fig. 1).
|
1–3 h by preventing free catalytic subunit from reassociating with the regulatory subunit. This intermediate-term prolongation of PKA activity contrasts with the long-term activation of PKA by Ap-Uch, which lasts for >20 h.
In the model, long-term PKA activity is assumed to correlate with LTM. This assumption appears plausible for the following reasons. Inhibition of PKA after 5-HT exposure blocks LTF, which is a correlate of LTM (28). Also, Ap-Uch mediates long-term PKA activation (28), and inhibition of Ap-Uch function or expression blocks LTF (25).
Methods of predicting nonlinear features of the model
The 15 differential equations (Eqs. 1–15, Table 1) can be represented as a vector equation,
 | (16) |
.
In Eq. 16, 
is a 15-dimensional vector of nonlinear functions; x = ([cAMP], [R], [C], [RC], [REG], [pREG], [mRNAREG], [Raf], [MAPKK], [MAPKKpp], [ERK], [ERKpp], Ppka, Perk, [Ap-Uch]) 
is the 15-dimensional vector of system variables; 
is the bifurcation parameter; 
denotes the time delay in Eqs. 2–4; and 
is the 38-dimensional vector of system parameters other than (i.e., the kinetic rate constants listed in Table 1). The full parameter space of the model, 
, is defined to include all real values for all parameters including . However, a restricted parameter space 
is also defined, which does not include , because a principle goal of bifurcation analysis is to subdivide this restricted parameter space into regions that exhibit qualitatively different dynamics when is varied. In practice, biochemical constraints (e.g., plausible rate constants and concentrations) constrain the portion of parameter space that is physiologically relevant. The concentration of 5-HT, [5-HT], is selected as the bifurcation parameter because [5-HT] is an important extracellular stimulus commonly varied in experiments, including experiments to determine a threshold for persistent PKA activation.
To use bifurcation and singularity theory to analyze the nonlinear properties of this dynamic system, steady states are first obtained by setting dx/dt = 0 in Eq. 16, yielding 15 nonlinear algebraic equations
 | (17) |
(in Eq. 16). Although has a significant effect on the transient dynamics of the system and might affect the stability of some steady states, it has no effect on the existence of a steady-state solution, or on the existence and location of limit points and singular points in the steady-state bifurcation diagrams. Therefore, setting the time delay to zero does not alter the steady-state bifurcation and singularity analysis described in this article. Hence, we simplify the steady-state bifurcation analysis by eliminating . We subsequently verify the predictions of the bifurcation analysis by simulation of the full model including the time delay.
Bifurcation diagrams, singularity theory, and classification diagrams
Two bifurcation diagrams are defined to be similar if the number, order, and orientation of the steady-state solutions x change in an identical way as the bifurcation parameter is varied. Singularity theory (8,29,30) was used to determine the regions in the restricted parameter space (not including ) that contain similar bifurcation diagrams, by tracing the region boundaries. Golubitsky and Schaeffer (8) showed that the boundaries consist of one of three types of hypersurfaces: a hysteresis surface (HS), an isola surface, or a double limit set (DLS). When the values of parameters are biochemically constrained to a certain range (e.g., nonnegative concentrations), a fourth type of hypersurface, a boundary limit set (BLS), subdivides the parameter space into allowed and disallowed regions. The hypersurfaces are defined within the full parameter space including (Eqs. 19–22, Appendix). However, their projections into the restricted parameter space are used to subdivide the restricted space into regions with qualitatively different bifurcation diagrams. DLS surfaces only exist with bifurcation diagrams displaying five or more steady-state solutions for specific values of the bifurcation parameter. The coexistence of more than three steady states was not observed in this study, and thus, DLS surfaces do not appear to exist in the parameter space of the present model.
Steady-state singular points or singularities are points in the full parameter space at which the Jacobian of the function f in Eq. 17 (see Appendix) has a zero eigenvalue (7,8). Singular points are defined to be of codimension n, if n–1 parameters, in addition to the bifurcation parameter, are specified to locate the point (see Appendix). A limit point, or saddle-node bifurcation point, is a point on a bifurcation diagram where the number of steady-state solutions changes by two as the bifurcation parameter is varied. A limit point is a singular point of codimension-1 because only the bifurcation parameter needs to be specified to determine the point (Appendix, Eq. 18). HS, isola, and DLS hypersurfaces are also comprised of singular points. These singular points are of codimension-2, because, to determine the point, two parameters ( and an additional parameter P1) are specified (Appendix, Eqs. 19, Eqs. 21, and Eqs. 22). Column 3 of Table 2 illustrates the types of bifurcation diagrams that exist at critical singular points on the various boundaries (HS, isola surface, or BLS) between the regions of the restricted parameter space. Columns 2 and 4 of Table 2 display the bifurcation diagrams obtained by universal unfolding of these hypersurfaces, with universal unfolding defined as perturbation away from the surface by a small change of a parameter other than (8,31). Unfolding HS leads to two possible scenarios, either a monotonically dependent bifurcation diagram with no limit point, or a hysteretic loop type bifurcation diagram (bistability) with two limit points that coalesce on the critical hysteresis surface (HS). Unfolding an isola surface causes an isolated branch of solutions to appear (see Table 2, Isola1 and Isola2). Table 2 illustrates two types of isolated branches, or isolae, i.e., elliptic and hyperbolic. Unfolding the BLS causes the limit point of the bifurcation diagram at the feasibility boundary of the bifurcation parameter to move either into or away from the feasibility range. All these unfolding scenarios are observed with the present model, as discussed below (Results, Fig. 3).
|
|
Computational methods
An outline of the general procedure for constructing a classification diagram is as follows (see Appendix for further details). A steady-state bifurcation diagram is constructed to identify codimension-1 limit points (e.g., Fig. 2 A). These limit points are used as initial points to construct a curve, or locus, of limit points. This locus constitutes a codimension-1 diagram (e.g., Fig. 2 B). An arc-length continuation scheme (32,33) is used to construct the curve (locus of limit points) by linear extrapolation of two points on the curve to determine a third point by the Newton-Raphson iteration scheme for solving nonlinear algebraic equations (34). The codimension-1 locus contains critical points that are of type codimension-2. The codimension-2 points are then used as initial points for constructing a diagram of loci of codimension-2 singularities, termed a classification diagram (e.g., Fig. 4). To construct this diagram, the arc-length continuation scheme is applied again.
|
|
0, as is illustrated in Figs. 5–8. The analysis, carried out on both positive and negative values of the bifurcation parameter before simplification to physiological values, is necessary to locate all steady-state solutions.
|
|
|
|
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
To further examine the steady-state solutions of the model, the PKA/Ap-Uch positive feedback loop (see Methods) was analyzed in greater detail. The synthesis of Ap-Uch plays a key role in regulating the strength of the positive feedback loop, because Ap-Uch is induced by PKA activation and subsequently hydrolyzes the R subunits of PKA to support long-term PKA activation. Using the two limit points (marked by crosses) in Fig. 2 A as initial codimension-1 points, a codimension-1 diagram was constructed in the parameter plane of the Ap-Uch synthesis rate constant (kApSyn) versus the bifurcation parameter [5-HT] (Fig. 2 B). To construct this diagram, the locus of limit points, also called the limit set (LS), was traced using an arc-length continuation scheme (see Methods). The two crosses marked on the LS in Fig. 2 B correspond to the two limit points in Fig. 2 A with kApSyn = 0.02 (standard value in Table 1). The LS in the (kApSyn,[5-HT]) plane is a closed loop.
Six codimension-2 singular points were found on LS, denoted in Fig. 2 B as HS1, HS2, Isola1, Isola2, BLS1, and BLS2. These singular points divide the range of kApSyn into seven intervals labeled 1–7 in Fig. 2 B. Each interval represents a distinct region of parameter space that is described by a different bifurcation diagram (see Fig. 3).
Bifurcation diagrams in each of these seven intervals were characterized. In interval 1 of Fig. 2 B where kApSyn is less than the critical value of HS1, PKA activity is a monotonic function of [5-HT] (bifurcation diagram in panel 1 of Fig. 3), because no limit point exists in this interval. PKA activity is at a very low level at [5-HT] = 0. In interval 2 in which kApSyn is greater than HS1, multiplicity appears. The bifurcation diagram exactly at the HS1 point is schematized in Table 2 (HS). After crossing HS1, the bifurcation diagram exhibits two limit points, a hysteresis loop, and an interval of [5-HT] with three steady states of PKA activity. A diagram from interval 2 is shown in panel 2 of Fig. 3, in which a hysteresis loop is illustrated by the transitions between the two branches of stable steady-state solutions, which occur upon reaching the two limit points.
Upon increasing kApSyn to intervals 3 and 4, another hysteresis surface point, HS2, is crossed. Two new limit points, and a narrow range of multiplicity between them, appear at negative [5-HT] (panels 3 and 4 of Fig. 3, insets). The bifurcation diagrams in panels 3 and 4 of Fig. 3, with four coexisting limit points, are of the mushroom type (8). In interval 4 (see Fig. 2 B, inset and the segment of the LS between BLS2 and HS1), a line with fixed kApSyn intersects the LS curve at four points, only one of which is at positive [5-HT]. In the corresponding bifurcation diagram (panel 4 of Fig. 3), only one limit point is at positive [5-HT]. By further increasing kApSyn in Fig. 2 B, HS1 and HS2 merge at the point (d: Isola1). Here, an isola hypersurface intersects the (kApSyn,[5-HT]) plane. The bifurcation diagram at Isola1 is illustrated in Table 2 (Isola1, hyperbolic), in which the two branches of the mushroom merge at the value of [5-HT] at Isola1. The number of limit points decreases from four to two. Upon crossing Isola1, kApSyn goes to intervals 5 and 6. An isola of steady-state solutions appears (panels 5 and 6 of Fig. 3). For those values of [5-HT] at which three steady states of PKA activity coexist, the lower and upper states are stable whereas the middle state is unstable (5–8). However, with an isola, the two stable steady-state branches are not connected to each other. Therefore, the system cannot pass from one branch to the other through a hysteresis mechanism or loop. However, isolae demonstrate a mechanism of irreversible, nonhysteretic transition between steady states. For example, suppose in either panel 5 or 6 of Fig. 3 that the system is initialized in the lower steady state. Then, upon increasing [5-HT], the lower steady state will disappear (right boundary of isola) and the system will transit to the upper steady state. The transition is irreversible in that if [5-HT] is then decreased to zero, the upper steady state remains stable, and PKA activity remains high. However, a brief but large dynamic perturbation of system variables or parameters might cause a state transition in either direction by forcing the system out of the basin of attraction of either steady state.
With a further increase of kApSyn, the isola shrinks to a point at Isola2 (Fig. 2 B, corresponding bifurcation diagram in Table 2 (Isola2, elliptic)). By crossing Isola2 to interval 7, the isola disappears, leaving only a unique steady state of high PKA activity (panel 7 of Fig. 3).
In addition to the effects of the HS and isola points (Fig. 2 B) on the model behavior, the effects of the boundary limit set (BLS) points also need to be considered to characterize the complete set of steady-state behaviors of the model. The BLS is the hypersurface in the full parameter space of limit points located at zero [5-HT] (i.e., at the physically meaningful boundary of [5-HT]). Fig. 2 B shows two codimension-2 singular points on the LS loop, denoted BLS1 and BLS2, at which the BLS hypersurface intersects the (kApSyn,[5-HT]) plane. At BLS1, one limit point in the 
bifurcation diagram is at zero [5-HT] (see BLS in Table 2). BLS1 divides the kApsyn interval containing a mushroom-type bifurcation diagram into two subdomains, intervals 3 and 4 in Fig. 2 B. Typical diagrams in these intervals are illustrated in panels 3 and 4 of Fig. 3, respectively. Upon unfolding BLS1 (i.e., perturbing kApSyn slightly above or below BLS1), either one steady-state solution (panel 3) or three solutions (panel 4) appear at [5-HT] = 0. At BLS2, the right limit point of the isola is at zero [5-HT], while the other limit point lies in negative [5-HT] (not shown). In Fig. 2 B, BLS2 divides the isola interval into two subdomains, intervals 5 and 6. Typical bifurcation diagrams are shown in panels 5 and 6 of Fig. 3, respectively. In Fig. 2 B, the curve that connects BLS2 and HS1 represents a threshold of 5-HT stimulation. If [5-HT] remains to the right of this curve, then the lower steady state of PKA activity no longer exists (the system is to the right of the lower limit points in Fig. 2 A or panels 2–5 of Fig. 3). The system will transit to the only stable state, of high PKA activity.
It is evident from the above that the LS in Fig. 2 B plays a key role in bifurcation analysis. The LS enables rapid determination of all qualitatively different bifurcation diagrams of PKA activity versus [5-HT] that can be obtained at different values of a system parameter (kApSyn). Furthermore, the LS allows for the identification of codimension-2 singular points (HS, isola, BLS points in Fig. 2 B). These points, and the associated hypersurfaces, are needed to characterize the effect of changing two or more parameters on the dynamics of the system. Specifically, loci of these codimension-2 points in the plane of any two parameters P1 and P2 (in the parameter space P) can be constructed. A classification diagram of these loci subdivides the plane into regions, each of which manifests a qualitatively different type of (x,) bifurcation diagram.
The loci of six codimension-2 singular points partition a key classification diagram into seven regions
In addition to kApSyn, a second parameter that determines the strength of the PKA/Ap-Uch positive feedback loop is kAp-Uch. The rate constant kAp-Uch describes degradation of the R subunit of PKA by Ap-Uch, enhancing the catalytic (C) activity of PKA (Fig. 1, Eqs. 2–4 in Table 1). Therefore, a classification diagram in the (kApSyn,kAp-Uch) plane, illustrating qualitative changes in dynamics as these parameters are varied, is likely to provide significant insight into regulation of the model's behavior. By continuation (varying kAp-Uch, see Methods) of the six codimension-2 singular points (in Fig. 2 B) in the (kApSyn,kAp-Uch) plane, this classification diagram of Fig. 4 was constructed. At any point on a codimension-2 curve, the bifurcation diagram of PKA activity versus [5-HT] is qualitatively like that of the corresponding codimension-2 singular point in Table 2. For example, at any point on the Isola1 curve, the bifurcation diagram will show an isola pinching off from another branch of stable solutions (Table 2). At any point within a region between curves in Fig. 4, the bifurcation diagram is qualitatively like the corresponding diagram in Fig. 3. For example, throughout region 2 of Fig. 4, the bifurcation diagram is like that of panel 2 of Fig. 3, with two limit points and a hysteresis interval of [5-HT] between the limit points. Three steady states will be present within this [5-HT] interval, and only one steady state will be present for all other values of [5-HT]. The specific diagrams of Fig. 3 lie along the line kAp-Uch = 0.007 min–1 in Fig. 4.
Only a subset of bifurcation diagrams is present for physiologically relevant 5-HT concentrations
By using the loci of codimension-1 singular points (LS in Fig. 2 B), different bifurcation diagrams (Fig. 3) can be derived and insights can be obtained into the ways in which these diagrams transit between each other via codimension-2 singular points (Table 2). Loci of these codimension-2 points partition classification diagrams (Fig. 4) into regions with distinct bifurcation diagrams. Analogous diagrams exist for any two parameters other than [5-HT]). These results were derived allowing [5-HT] to vary over negative, zero, and positive values. Physiologically, only nonnegative [5-HT] is relevant. Can the portions of bifurcation diagrams at negative [5-HT] be neglected and a simplified analysis derived that focuses on physiological significance?
To address this question, the seven bifurcation diagrams of Fig. 3 were further inspected. By neglecting the portions of diagrams with negative [5-HT], the diagrams can be reclassified into four types. In panel 1 of Fig. 3, the PKA activity exhibits a unique, monotonically increasing function for [5-HT] 
, but because the PKA/Ap-Uch positive feedback is weak (low kApSyn), PKA activity remains low. In panels 2 and 3 of Fig. 3, PKA activity manifests three steady-state values within a positive range of [5-HT]. In panels 4 and 5 of Fig. 3, PKA activity manifests three steady states within a range that includes [5-HT] = 0. In panels 6 and 7 of Fig. 3, a unique steady state again exists, but because the PKA/Ap-Uch positive feedback is strong (high kApSyn), PKA activity is high even at zero [5-HT]. These bifurcation diagrams with [5-HT] are plotted in Fig. 5. Fig. 5, A–D, corresponds to panels 1–7 of Fig. 3 as follows: Fig. 5 A to panel 1; Fig. 5 B to panels 2 and 3; Fig. 5 C to panels 4 and 5; and Fig. 5 D to panels 6 and 7.
There are three codimension-2 singular points with either positive or zero [5-HT] in Fig. 2 B: HS1, BLS1, and BLS2. These three points divide the range of kApSyn into four intervals, a, b, c, and d. The intervals a–d correspond to the intervals 1–7 of Fig. 2 B as follows: a to 1; b to 2 and 3; c to 4 and 5; and d to 6 and 7. The corresponding bifurcation diagram in interval a, below HS1, is shown in Fig. 5 A. Interval b, bounded by HS1 and BLS1, has the bifurcation diagram shown in Fig. 5 B. In interval b, a line with fixed kApSyn will intersect the LS in Fig. 2 B at two points for [5-HT] 
0. These points correspond to the two limit points in the bifurcation diagram shown in Fig. 5 B (LP1 and LP2). In interval c of kApSyn, bounded by BLS1 and BLS2, the bifurcation diagram corresponds to Fig. 5 C. A line with fixed kApSyn in Fig. 2 B intersects the LS at only one point for [5-HT] 0, and the corresponding bifurcation diagram in Fig. 5 C also shows one limit point. In interval d, with kApSyn above BLS2, the bifurcation diagram corresponds to Fig. 5 D and has no limit point for [5-HT] 0.
The bifurcation diagram in Fig. 5 B illustrates a hysteretic loop between two limit points. Such a structure is termed a reversible switch or toggle switch, because alterations in [5-HT] can carry the system to the left or right of both limit points, inducing transitions of PKA activity between the upper and lower steady states. This multiplicity is termed bistability since the upper and lower solutions are stable whereas the middle solution is unstable. The limit point with larger [5-HT], LP1, defines a threshold of [5-HT] required to induce high PKA activity. Physiologically, such a threshold might act as a filter, by which a weak 5-HT stimulus can be regarded as noise and neglected. A saturation of PKA activity was also observed in this bifurcation diagram for [5-HT] > 0.16. Saturation might help to protect neurons against excessive gene induction and energy expenditure in response to frequent stimuli.
Although Fig. 5 loses part of the detailed information about how these bifurcation diagrams originate (Fig. 3) and transit between each other, it is derived from the overall bifurcation and singularity theory, and focuses on the physiologically relevant concentrations of 5-HT. A classification diagram similar to Fig. 4, but restricted to [5-HT] 0, can be constructed (Fig. 6) by continuing the codimension-2 singular points of Fig. 2 B that have nonnegative [5-HT] (HS1, BLS1, and BLS2) in the kApSyn-kAp-Uch plane. Regions a–d bounded by the loci of these singular points have the bifurcation diagrams illustrated in Fig. 5, A–D, respectively.
|
In the classification diagram of Fig. 6, region c, which corresponds to the irreversible-switch bifurcation diagram, is relatively large. This observation suggests that, generically, the irreversible switch will be preserved for small parameter variations. However, the (kApSyn,kAp-Uch) plane is only one of the possible slices through the restricted parameter space. Does the irreversible-switch region retain substantial size if the parameter space is examined using different parameters? To address this question, the (kphos1,kphos2) parameter plane was used to construct an analogous classification diagram (Fig. 7) (kphos1 is the rate constant with which PKA phosphorylates TF-1; kphos2 is the rate constant with which ERK phosphorylates TF-2). Although Figs. 6 and 7 appear somewhat different, their essential characteristics are similar. In both figures, one hysteretic surface (HS1) and two boundary limit sets (BLS1 and BLS2) divide the plane into four regions that correspond to the topologically distinct bifurcation diagrams in Fig. 5, A–D. As in Fig. 6, the region corresponding to the irreversible switch (region c in Fig. 7) remains relatively large. The other regions a, b, and d also remain of substantial size. This result strengthens the conclusion that the irreversible switch, as well as the dynamics represented in the other regions, are robust to small parameter variations.
To further assess the robustness of the irreversible switch, the effects of some large parameter changes on region c, Figs. 6 and 7, were examined. Fig. 8 A illustrates the consequences of doubling either kphos1 or kphos2. It is evident that the size of region c decreases somewhat upon doubling kphos1 or kphos2. The shift of the BLS2 boundary to the left is more significant than that of BLS1. However, the region is still substantially large. In particular, when kAp-Uch is high, a large range of kApSyn generates an irreversible switch. Fig. 8 B exhibits the consequences of doubling either kAp-Uch or kApSyn. Increasing kApSyn does not significantly alter the size of the irreversible switch region bounded by BLS1 and BLS2. Increasing kAp-Uch can considerably enlarge the irreversible switch domain. These results illustrate that the irreversible switch exists in a substantial region of parameter space, and is robust to large parameter changes. This robustness suggests that the irreversible switch mechanism might be observed to function in vivo.
Simulations verify steady-state bifurcation analysis
The above bifurcation analyses do not consider the time delay (Eq. 16) in the original model of Fig. 1. The presence of a time delay cannot eliminate any of the steady states in Figs. 2–8
. However, a delay can affect the dynamic behavior as the system evolves toward steady states, and may alter their stability. Therefore, it is important to complement the steady-state analyses with simulations in which the system is initially far from steady state.
Four distinct steady-state bifurcation diagrams are plotted in Fig. 9. As kApSyn is increased, a monotonic steady state of low PKA activity (kApSyn = 0.0022) changes to a hysteretic loop with two limit points and three steady states of PKA activity between these points (kApSyn = 0.021). At higher kApSyn (0.03), an irreversible switch is obtained, with only one limit point at positive [5-HT] and three steady states at [5-HT] = 0. At yet higher kApSyn (0.1), only one steady state exists, with high PKA activity. These four bifurcation curves correspond to the four regions of Fig. 6 and are marked in Fig. 6 by the dots on the line kAp-Uch = 0.007 min–1.
Fig. 9, B–E, illustrate simulations of the model (Table 1) with the time delay ( = 250 min). By applying a long (100 h) stimulus, these simulations verify the stability of the steady states in Fig. 9 A. With kApSyn = 0.0022, 100 h of 10 µM 5-HT stimulation elevates PKA activity to a plateau (Fig. 9 B). The activity returns to its baseline rapidly after the stimulus ends, and this baseline corresponds to the lowest steady state in Fig. 9 A with [5-HT] = 0. With kApSyn = 0.021 (the reversible switch bifurcation diagram in Fig. 9 A), similar dynamics to Fig. 9 B are observed (Fig. 9 C). However, after the same 5-HT stimulation, PKA activity slowly returns to the steady-state value (in Fig. 9 A with [5-HT] = 0, kApSyn = 0.021). With kApSyn = 0.03 (the irreversible switch bifurcation diagram in Fig. 9 A), a prolonged stimulus of 10 µM 5-HT induces a transition to a high, stable level of PKA activity. This state remains stable when [5-HT] returns to 0 µM at t = 100 h, and PKA activity is permanently locked in a high state (Fig. 9 D). Finally, with very strong positive feedback (kApSyn = 0.1) (Fig. 9 E), PKA activity converges to the highest stable state of Fig. 9 A even with [5-HT] fixed at 0 µM.
The model (Table 1) contains three delay differential equations (DDEs) (Eqs. 2–4). To locate the steady-state solutions, limit points, and singular points, the model in vector form (Eq. 16) was reduced to the algebraic system (Eq. 17) without time delay. Time delays do not alter the existence of a steady state, or the locations of codimension-1 and -2 steady-state singular points (35). The focus of this article is to use steady-state singularity analysis to compute the singular points and construct classification diagrams of parameter space to determine the regions in which different steady-state bifurcation diagrams exist. This analysis is invariant to the presence of a time delay.
Analyses of DDEs and ODEs do differ in terms of the procedure to determine the stability of a steady state. Time delays may destabilize a steady state under certain conditions. In this study, we did not systematically investigate how the time delay affects the stability of steady states through analysis of the eigenvalues of linearized DDEs. Instead, we verified the stability of several key steady-state solutions on the different bifurcation diagrams in Fig. 9 A by simulation of the model using different values of . For example, we investigated the stability of the upper steady-state branch of the irreversible switch at [5-HT] = 0 (with kApSyn = 0.03) in Fig. 9 A. A delay equal to or smaller than that used by Pettigrew et al. (17) ( 
250 min) has a negligible impact on the transients as the system converges to the stable steady state. Increasing the delay 20-fold to 5000 min causes the system to converge to the steady state by way of damped oscillations with small amplitude. With yet a longer time delay (10,000 min), the amplitude of the damped oscillation increases, but the stability is still not altered. Similar results were found for the lower steady-state branch of the irreversible switch (with kApSyn = 0.03) in Fig. 9 A. The middle steady-state branch is unstable, therefore it cannot be reached with transient simulations. We also examined the stability of steady-state solutions of the three other bifurcation diagrams in Fig. 9 A and found that the stability is not altered by time delays as large as 5000 min. Therefore, we conclude that physiologically reasonable time delays do not alter the stability of the steady states in the parameter space that we investigated.
By applying pulsed stimulation protocols, which are often implemented experimentally (27,28), simulations were performed in the parameter regimes of the reversible switch (kApSyn = 0.021) and irreversible switch (kApSyn = 0.03) (Fig. 10). After five pulses of 5-HT (10 µM for 5 min, with an interstimulus interval of 20 min), PKA activity is strongly activated at 24 h in both cases (Fig. 10 A). However, PKA activity gradually drops to its baseline in the reversible switch case after removing 5-HT, whereas PKA activity remains at high level in the irreversible switch case (Fig. 10 B). After three pulses of 5-HT, PKA activity is only slightly and briefly elevated in both cases (not shown). Thus, three pulses of 5-HT are subthreshold, whereas five pulses cross the threshold for long-term PKA activation. Empirically, five pulses of 5-HT, but not three, induce long-term PKA activation (36), as well as LTF (10,37). In vivo, four or five spaced shocks to Aplysia induce long-term sensitization (LTS) of withdrawal reflexes, whereas three shocks do not (11). These observations suggest the threshold stimulus duration for long-term PKA activation (five pulses) corresponds to a threshold for the formation of LTS, a simple form of LTM.
|
7.5 µM [5-HT]. This stimulus is the minimum amplitude to induce high PKA activity by using the five-pulse stimulation protocol. However, the steady-state threshold for PKA activity in Fig. 9 A is just above a 5-HT stimulus of 0.106 µM. This steady-state threshold is the minimal stimulus amplitude to induce high PKA activity by continuous 5-HT application or very long stimulus duration. A simulation with [5-HT] held at 0.15 µM, slightly above the steady-state threshold, showed that to induce high PKA activity, 4500 min of stimulation was required. |
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONAPPENDIXACKNOWLEDGEMENTSREFERENCES |
|---|
0 in Fig. 2 B. It might appear that if analysis focuses on physiologically relevant parameter values, a subset of the singular points may suffice to classify parameter space. However, for a complete analysis it is necessary to investigate both the positive and negative ranges of the bifurcation parameter [5-HT]. This analysis is based on locating the limit points of the model, some of which are at negative [5-HT] (Fig. 2 B), to characterize the possible bifurcation diagrams. The limit points are used to locate the codimension-2 singular points (Fig. 2 B), which mark the boundaries between various types of bifurcation diagrams (Fig. 3). In our analysis, most of these singular points are at negative [5-HT] (Fig. 2 B). Only after identifying these points, and the intervals between them, can we classify the different types of bifurcation diagrams and construct specific examples (Fig. 3). If only positive values of the bifurcation parameter are examined, we risk missing a boundary between regions in the full classification diagram (Fig. 4), as well as incorrectly identifying the transitions between them (see Table 2). For example, examining negative [5-HT] allows us to locate the BLS2 limit point. BLS2 represents the boundary between regions c and d in the classification diagram (Fig. 6) corresponding to the bifurcation diagrams of the irreversible switch type (Fig. 5 C) and monotonic type (Fig. 5 D), respectively. The irreversible switch is a potential mechanism for the induction or maintenance of LTM. In addition, without the analysis with negative [5-HT], it is not possible to establish that the bifurcation diagrams containing an irreversible switch (Fig. 5 C, corresponding to both intervals 4 and 5 of kApSyn in Fig. 2 B) mainly originate from a broken isola bifurcation, and not from a broken hysteretic loop. The origin of the irreversible switch can be observed from Figs. 2 B and 4, in which interval or region 5 (corresponding to the isola bifurcation diagram in Fig. 3, panel 5) is large, whereas interval or region 4 (corresponding to the hysteretic bifurcation diagram in Fig. 3, panel 4) is small (see Fig. 2 B, inset).
Bifurcation analysis and simulation complement each other
Simulation is essential to study the temporal evolution of a model, whereas bifurcation analysis provides a steady-state view. These two methods of analysis complement each other, as illustrated by comparing the simulations of Figs. 9, B–E, and 10 with the bifurcation diagrams of Fig. 9 A. Simulations with the full model (Table 1) including the time delay are consistent with steady-state features of bifurcation diagrams derived without . This comparison emphasizes the necessity of performing both steady-state analysis and simulations to understand the complex interactions within a biochemical network. Bifurcation analysis systematically characterizes the parameter space to determine regions that generate qualitatively different steady-state or oscillatory solutions. For complex models, this analysis is necessary, because it is computationally impractical to characterize the dynamics by a very large number of random or systematic simulations throughout the multidimensional parameter space. In contrast, simulations can generate time courses of system variables that provide detailed information on how the variables respond to stimuli and interact with each other. Simulations also provide insights into the methods that may be used to reach steady states (e.g., by determining the minimum stimulus amplitude and duration required to reach a specific steady state).
Irreversible switch dynamics may contribute to the formation of LTM
Irreversible switches have been suggested to be important for cellular differentiation and development (18–23); entry into mitosis in the cell cycle (3); ligand-induced autocrine growth factor signaling (38); Sonic Hedgehog signaling (39); and long-term memory formation (40–44). The common feature is conversion of a transient stimulus into an irreversible response through the activation of a network of intrinsically reversible signaling molecules and gene products.
Several types of irreversible switches have been proposed to contribute to the maintenance of LTM. Lisman and Fallon (41) and Roberson and Sweatt (44,45) suggested several possible switches that might sustain long-term changes in the strength of synapses. Lisman and co-workers (43,46,47) proposed CaMKII activity may be sustained indefinitely by positive feedback in which active, phosphorylated CAMKII multimers further phosphorylate their subunits, and this CAMKII activity might maintain long-lasting synaptic strengthening or potentiation (LTP) and therefore LTM. Iyengar and co-workers (40,48) showed that a feedback loop involving protein kinase activities and second messengers (including protein kinase C, Raf, MAP kinase kinase, and MAP kinase) can produce a bistable molecular switch at a synapse, which when activated may induce LTP. Kandel and co-workers proposed a positive feedback mechanism for protein oligomerization (49–51) that may be important for sustaining persistent modification of synaptic structure, allowing for LTM. The irreversible switch exhibited by the present model is similarly based on a positive feedback loop, with reciprocal induction of PKA activation and Ap-uch expression (see Methods). Such a switch (Figs. 5 C and 6) could maintain long-lasting PKA activation, transcription factor (TF) phosphorylation, and gene induction in response to relatively brief stimuli, and could therefore be important for consolidation of LTM, which is known to require gene expression.
Switches based on positive feedback loops, such as the irreversible switch (Fig. 5 C), are commonly described by deterministic ODEs or DDEs, which use chemical concentrations as continuous variables. This representation implicitly assumes that the number of molecules in the system is sufficiently large to ensure that deterministic, mass-action kinetics can describe the dynamics. However, if the copy numbers of some key kinases, or other molecular species, in the feedback loop are relatively low (e.g., tens), then stochastic noise and fluctuations in molecular copy numbers may induce spontaneous transitions between states (23,52). Therefore, molecular noise may be a key factor that affects the robustness of the bistable switch when the copy numbers of molecules are low.
Whether the PKA/Ap-Uch positive feedback loop does, in fact, play an important role in the formation of LTM remains to be determined experimentally. If this loop exhibits switchlike behavior, then a relatively brief activation of PKA by an agonist that elevates cAMP levels, such as forskolin or 5-HT, should induce a subsequent, long-lasting (many hours or days) phase of PKA activation and Ap-uch expression. Indeed, in Aplysia neuronal ganglia, application of five 5-min pulses of 5-HT induces PKA activation measured 24 h later (36). Although this result is compatible with long-lasting positive feedback, the time course of PKA activation needs to be examined at additional time points. Crucially, concomitant long-lasting induction of Ap-uch would need to occur to support the hypothesis that brief stimuli can induce irreversible switchlike activation of the positive feedback. Furthermore, a molecular pathway from PKA activation to Ap-uch gene induction should be identified. As discussed above (see Methods), Ap-uch does not appear to be induced via phosphorylation of the CREB1 transcriptional activator. However, PKA activation could lead to phosphorylation of an unspecified TF that induces Ap-uch. A major goal of this article, not dependent on these model details, is to point out that positive feedback loops may be important for sustaining long-lasting kinase activation and gene expression to induce, or maintain, LTM.
We note that in vivo, any biochemical switch associated with the induction of LTM is not likely to remain irreversible. It is likely that processes of molecular turnover, not represented in the present model, would tend to return kinase activities and synaptic strengths to average levels. For example, several current theories of long-term synaptic plasticity envision homeostatic mechanisms that rescale synaptic weights over time (53,54).
Classification diagrams can help predict effects of drugs or other environmental stimuli
Classification diagrams of parameter space (e.g., Figs. 4 and 6–8) provide boundaries within the parameter space that differentiate among multiple bifurcation diagrams and multiple stimulus-response outcomes (e.g., hysteresis or irreversible switches). Classification diagrams help to determine which dynamic features are robust to variation of parameter values. Only those bifurcation diagrams that occupy a region of significant size in multiple classification diagrams (e.g., Figs. 6 c and 7 c, and perturbations of this region in Fig. 8) are likely to describe the physiological dynamics.
Codimension-1 bifurcation diagrams analogous to Fig. 2 B can be used to determine parameters that, if varied by pharmacological or environmental intervention, would have a major qualitative effect on the dynamics of the system. For example, Fig. 2 B illustrates that a compound affecting the rate of Ap-uch transcription or protein synthesis (the rate constant kApSyn) would significantly alter the dynamics (e.g., by changing hysteretic behavior to monotonic or irreversible-switch behavior). Diagrams similar to Fig. 2 B can be combined with considerations of which parameters could be altered relatively selectively by chemical or other stimuli. In this way, potential targets could be identified for drug design and therapy, as well as for toxins or other environmental stimuli.
Classification diagrams analogous to Figs. 4 and 6 provide a visual map to determine how the system will respond if two parameters are varied simultaneously by pharmacological or environmental stimuli, or if one parameter is varied in systems characterized by differing values of a second parameter. Therefore, classification diagrams may provide a visual means of determining the effect of two drugs applied simultaneously, and possibly illustrate a synergetic effect, which would not be observed if the drugs were applied individually or sequentially. For example, Fig. 6 illustrates that compounds or drugs might enhance the induction of LTM by altering either kApSyn or kAp-Uch or both. Increasing the strength of positive feedback by increasing kApSyn and/or kAp-Uch favors more prolonged and greater activation of PKA. A monotonic relationship of PKA activity versus [5-HT] (region a) is converted sequentially to a hysteretic reversible switch (region b), then to an irreversible switch (region c). In vivo, compounds that increase Ap-Uch synthesis (kApSyn) and/or increase the effectiveness with which Ap-Uch promotes PKA activation (kAp-Uch) would plausibly promote persistent PKA activation, TF phosphorylation, and synaptic strengthening. Near boundary curves of classification diagrams, the system is particularly sensitive to parameter changes, because small changes can qualitatively alter system behavior.
Pettigrew et al. (17) proposed (their Fig. 6 A) that a compound able to reduce the dephosphorylation rates of MAPKK (kb, MAPKK) and ERK (kb, ERK), might enhance 5-HT-induced Ap-uch expression and PKA activation. The effect of reducing both these rate constants was therefore examined. Bifurcation analysis of the unaltered system (kb, MAPKK = kb, ERK = 0.12 min–1) and the altered system (kb, MAPKK = kb, ERK = 0.06 min–1) (other parameters as in Table 1) show that the qualitative nature of the bifurcation diagram (a hysteretic loop/reversible switch) is not changed by decreasing kb, MAPKK and kb, ERK, although the threshold for PKA activation by 5-HT is considerably reduced (not shown). An irreversible-switch type of bifurcation diagram cannot be generated by further decreasing kb,MAPKK and kb,ERK (not shown), because the strength of the positive feedback loop (which is governed by kApSyn and kAp-Uch) remains insufficient. These simulations suggest alteration of kinetic parameters outside the positive feedback loop can decrease the threshold for 5-HT necessary for significant PKA activation, but is less likely to enhance the duration of PKA activation, TF phosphorylation, and gene induction. The analysis predicts that alteration of kinetic parameters within the positive feedback loop might be more likely to enhance both the intensity and duration of gene induction, and the
