We describe new methods for initializing the computation of homoclinic orbits for maps in a state space with arbitrary dimension and for detecting their bifurcations. The initialization methods build on known and improved methods for computing one-dimensional stable and unstable manifolds. The methods are implemented in MatcontM, a freely available toolbox in Matlab for numerical analysis of bifurcations of ï¬xed points, periodic orbits and connecting orbits of smooth nonlinear maps. The bifurcation analysis of homoclinic connections under variation of one parameter is based on continuation methods and allows to detect all known codimension 1 and 2 bifurcations in 3D maps, including tangencies and generalized tangencies. MatcontM provides a graphical user interface, enabling interactive control for all computations. As the prime new feature we discuss an algorithm for initializing connecting orbits in the important special case where either the stable or unstable manifold is one-dimensional, allowing to compute all homoclinic orbits to saddle points in three-dimensional maps. We illustrate this algorithm in the study of the adaptive control map, a 3D map introduced in 1991 by Frouzakis, Adomaitis and Kevrekidis, to obtain a complete bifurcation diagram of the resonance horn in a 1:5 Neimark-Sacker bifurcation point, revealing new features.
The current opinion in epilepsy surgery is that successful surgery is about removing pathological cortex in the anatomic sense. This contrasts with recent developments in epilepsy research where epilepsy is seen as a network disease. Computational models offer a framework to investigate the influence of networks, as well as local tissue properties, and to explore alternative resection strategies. Here we study, using such a model, the influence of connections on seizures and how this might change our traditional views of epilepsy surgery. We use a simple network model consisting of four interconnected neuronal populations. One of these populations can be made hyperexcitable, modeling a pathological region of cortex. Using model simulations the effect of surgery on the seizure rate is studied. We find that removal of the hyperexcitable population is, in most cases, not the best approach to reduce the seizure rate. Removal of normal populations located at a crucial spot in the network, the ''driver'', is typically more effective in reducing seizure rate. This work strengthens the idea that network structure and connections may be more important than localizing the pathological node. This can explain why lesionectomy may not always be sufficient.
Objective: In postanoxic coma, EEG patterns indicate the severity of encephalopathy and typically evolve in time. We aim to improve the understanding of pathophysiological mechanisms underlying these EEG abnormalities.
Methods: We used a mean field model comprising excitatory and inhibitory neurons, local synaptic connections, and input from thalamic afferents. Anoxic damage is modeled as aggravated short-term synaptic depression, with gradual recovery over many hours. Additionally, excitatory neurotransmission is potentiated, scaling with the severity of anoxic encephalopathy. Simulations were compared with continuous EEG recordings of 155 comatose patients after cardiac arrest.
Results: The simulations agree well with six common categories of EEG rhythms in postanoxic encephalopathy, including typical transitions in time. Plausible results were only obtained if excitatory synapses were more severely affected by short-term synaptic depression than inhibitory synapses.
Conclusions: In postanoxic encephalopathy, the evolution of EEG patterns presumably results from gradual improvement of complete synaptic failure, where excitatory synapses are more severely affected than inhibitory synapses. The range of EEG patterns depends on the excitation-inhibition imbalance, probably resulting from long-term potentiation of excitatory neurotransmission.
Healthy or pathological states of nociceptive subsystems determine different stimulus-response relations measured from quantitative sensory testing. In turn, stimulus-response measurements may be used to assess these states. In a recently developed computational model, six model parameters characterize activation of nerve endings and spinal neurons. However, both model nonlinearity and limited information in yes-no detection responses to electrocutaneous stimuli challenge to estimate model parameters. Here, we address the question whether and how one can overcome these difficulties for reliable parameter estimation. First, we fit the computational model to experimental stimulus-response pairs by maximizing the likelihood. To evaluate the balance between model fit and complexity, i.e., the number of model parameters, we evaluate the Bayesian Information Criterion. We find that the computational model is better than a conventional logistic model regarding the balance. Second, our theoretical analysis suggests to vary the pulse width among applied stimuli as a necessary condition to prevent structural non-identifiability. In addition, the numerically implemented profile likelihood approach reveals structural and practical non-identifiability. Our model-based approach with integration of psychophysical measurements can be useful for a reliable assessment of states of the nociceptive system.
To continue a branch of homoclinic solutions starting from a Bogdanov-Takens (BT) point in parameter and state space, one needs a predictor based on asymptotics for the bifurcation parameter values and the corresponding small homoclinic orbits in the phase space. We derive two explicit asymptotics for the homoclinic orbits near a generic BT point. A recent generalization of the Lindstedt-Poincare (L-P) method is applied to approximate a homoclinic solution of a strongly nonlinear autonomous system that results from blowing up the BT normal form. This solution allows us to derive an accurate second-order homoclinic predictor to the homoclinic branch rooted at a generic BT point of an n-dimensional ordinary differential equation (ODE). We prove that the method leads to the same homoclinicity conditions as the classical Melnikov technique, the branching method, and the regular perturbation (R-P) method. However, it is known that the R-P method leads to a `parasitic turn' near the saddle point. The new asymptotics based on the L-P method do not have this turn, making them more suitable for numerical implementation. We show how to use these asymptotics to calculate the initial data to continue homoclinic orbits in two free parameters. The new homoclinic predictors are implemented in the MATLAB continuation package MatCont to initialize the continuation of homoclinic orbits from a BT point. Two examples with multidimensional state spaces are included.
The second-order predictor for the homoclinic orbit is applied to the Gray-Scott model. The problem is used to illustrate the approximation of the homoclinic orbits near a generic Bogdanov-Takens bifurcation in n-dimensional systems of differential equations. In the process we show that it is necessary to take (usually ignored) cubic terms in the Bogdanov-Takens normal form into account to derive a correct second-order prediction for the homoclinic bifurcation curve. The analytic solutions are compared with those obtained by numerical continuation.
Sensitization is an example of malfunctioning of the nociceptive pathway in either the peripheral or central nervous system. Using quantitative sensory testing, one can only infer sensitization, but not determine the defective subsystem. The states of the subsystems may be characterized using computational modeling together with experimental data. Here, we develop a neurophysiologically plausible model replicating experimental observations from a psychophysical human subject study. We study the effects of single temporal stimulus parameters on detection thresholds corresponding to a 0.5 detection probability. To model peripheral activation and central processing, we adapt a stochastic drift-diffusion model and a probabilistic hazard model to our experimental setting without reaction times. We retain six lumped parameters in both models characterizing peripheral and central mechanisms. Both models have similar psychophysical functions, but the hazard model is computationally more efficient. The model-based effects of temporal stimulus parameters on detection thresholds are consistent with those from human subject data.
Measurements of neuronal signals during human seizure activity and evoked epileptic activity in experimental models suggest that, in these pathological states, the individual nerve cells experience an activity driven depolarization block, i.e. they saturate. We examined the effect of such a saturation in the Wilson-Cowan formalism by adapting the nonlinear activation function; we substituted the commonly applied sigmoid for a Gaussian function. We discuss experimental recordings during a seizure that support this substitution. Next we perform a bifurcation analysis on the Wilson-Cowan model with a Gaussian activation function. The main effect is an additional stable equilibrium with high excitatory and low inhibitory activity. Analysis of coupled local networks then shows that such high activity can stay localized or spread. Specifically, in a spatial continuum we show a wavefront with inhibition leading followed by excitatory activity. We relate our model simulations to observations of spreading activity during seizures.
In this paper we develop an extended center manifold reduction method: a methodology to analyze the formation and bifurcations of small-amplitude patterns in certain classes of multi-component, singularly perturbed systems of partial differential equations. We specifically consider systems with a spatially homogeneous state whose stability spectrum partitions into eigenvalue groups with distinct asymptotic properties. One group of successive eigenvalues in the bifurcating group are widely interspaced, while the eigenvalues in the other are stable and cluster asymptotically close to the origin along the stable semi-axis. The classical center manifold reduction provides a rigorous framework to analyze destabilizations of the trivial state, as long as there is a spectral gap of sufficient width. When the bifurcating eigenvalue becomes commensurate to the stable eigenvalues clustering close to the origin, the center manifold reduction breaks down. Moreover, it cannot capture subsequent bifurcations of the bifurcating pattern. Through our methodology, we formally derive expressions for low-dimensional manifolds exponentially attracting the full flow for parameter combinations that go beyond those allowed for the (classical) center manifold reduction, i.e. to cases in which the spectral gap condition no longer can be satisfied. Our method also includes an explicit description of the flow on these manifolds and thus provides an analytical tool to study subsequent bifurcations. Our analysis centers around primary bifurcations of transcritical type -- that can be either of co-dimension 1 or 2 -- in two- and three-component PDE systems. We employ our method to study bifurcation scenarios of small-amplitude patterns and the possible appearance of low-dimensional spatio-temporal chaos. We also exemplify our analysis by a number of characteristic reaction-diffusion systems with disparate diffusivities.
Objective. Continuous application of high-frequency deep brain stimulation (DBS) often effectively reduces motor symptoms of Parkinson's disease patients. While there is a growing need for more effective and less traumatic stimulation, the exact mechanism of DBS is still unknown. Here, we present a methodology to exploit the plasticity of GABAergic synapses inside the external globus pallidus (GPe) for the optimization of DBS. Approach. Assuming the existence of spike-timing-dependent plasticity (STDP) at GABAergic GPe-GPe synapses, we simulate neural activity in a network model of the subthalamic nucleus and GPe. In particular, we test different DBS protocols in our model and quantify their influence on neural synchrony. Main results. In an exemplary set of biologically plausible model parameters, we show that STDP in the GPe has a direct influence on neural activity and especially the stability of firing patterns. STDP stabilizes both uncorrelated firing in the healthy state and correlated firing in the parkinsonian state. Alternative stimulation protocols such as coordinated reset stimulation can clearly profit from the stabilizing effect of STDP. These results are widely independent of the STDP learning rule. Significance. Once the model settings, e.g., connection architectures, have been described experimentally, our model can be adjusted and directly applied in the development of novel stimulation protocols. More efficient stimulation leads to both minimization of side effects and savings in battery power.
We analyze two ordinary differential equation (ODE) models for atherosclerosis. The ODE models describe long time evolution of plaques in arteries. We show how the dynamics of the first atherosclerosis model (model A) can be understood using codimension-two bifurcation analysis. The Low-Density Lipoprotein (LDL) intake parameter (d) is the first control parameter and the second control parameter is either taken to be the conversion rate of macrophages (b) or the wall shear stress (σ). Our analysis reveals that in both cases a Bogdanov-Takens (BT) point acts as an organizing center. The bifurcation diagrams are calculated partly analytically and to a large extent numerically using AUTO07 and MATCONT. The bifurcation curves show that the concentration of LDL in the plaque as well as the monocyte and the macrophage concentration exhibit oscillations for a certain range of values of the control parameters. Moreover, we find that there are threshold values for both the cholesterol intake rate dcritdcrit and the conversion rate of the macrophages bcritbcrit, which depend on the values of other parameters, above which the plaque volume increases with time. It is found that larger conversion rates of macrophages lower the threshold value of cholesterol intake and vice versa. We further argue that the dynamics for model A can still be discerned in the second model (model B) in which the slow evolution of the radius of the artery is coupled self-consistently to changes in the plaque volume. The very slow evolution of the radius of the artery compared to the other processes makes it possible to use a slow manifold approximation to study the dynamics in this case. We find that in this case the model predicts that the concentrations of the plaque constituents may go through a period of oscillations before the radius of the artery will start to decrease. These oscillations hence act as a precursor for the reduction of the artery radius by plaque growth.
An improved homoclinic predictor at a generic codim 2 Bogdanov-Takens (BT) bifucation is derived. We use the classical "blow-up" technique to reduce the canonical smooth normal form near a generic BT bifurcation to a perturbed Hamiltonian system. With a simple perturbation method, we not only rederive the known approximation for the bifurcation parameter values but also obtain an explicit first-order correction of the unperturbed homoclinic orbit. To obtain the normal form on the center manifold, we apply the standard parameter-dependent center manifold reduction combined with the normalization, that is based on the Fredholm solvability of the homological equation. However, by systematically solving all linear systems appearing from the homological equation, we remove an ambiguity in the parameter transformation existing in the literature. The actual implementation of the improved predictor in MatCont and numerical examples illustrating its effciency are discussed.
We consider travelling waves (fronts, pulses and periodics) in spatially extended one dimensional neural field models. We demonstrate for an excitatory field with linear adaptation that, in addition to an expected stable pulse solution, a stable anti-pulse can exist. Varying the adaptation strength we unravel the organizing centers of the bifurcation diagram for fronts and pulses, with a mixture of exact analysis for a Heaviside firing rate function and novel numerical schemes otherwise. These schemes, for non-local models with space-dependent delays, further allow for the construction and continuation of periodic waves. We use them to construct the dispersion curve -- wave speed as a function of period -- and find that they can be oscillatory and multi-valued, suggesting bistability of periodic waves. A kinematic theory predicts the onset of wave instabilities at stationary points in the dispersion curve, leading to period doubling behaviour, and is confirmed with direct numerical simulations. We end with a discussion of how the construction of dispersion curves may allow a useful classification scheme of neural field models for epileptic waves.
At one level of abstraction neural tissue can be regarded as a medium for turning local synaptic activity into output signals that propagate over large distances via axons to generate further synaptic activity that can cause reverberant activity in networks that possess a mixture of excitatory and inhibitory connections. This output is often taken to be a firing rate, and the mathematical form for the evolution equation of activity depends upon a spatial convolution of this rate with a fixed anatomical connectivity pattern. Such formulations often neglect the metabolic processes that would ultimately limit synaptic activity. Here we reinstate such a process, in the spirit of an original prescription by Wilson and Cowan [1972, ``Excitatory and inhibitory interactions in localised populations of model neurons'', Biophysical Journal 12, 1-24.], using a term that multiplies the usual spatial convolution with a moving time average of local activity over some refractory time-scale. This modulation can substantially affect network behaviour, and in particular give rise to periodic travelling waves in a purely excitatory network (with exponentially decaying anatomical connectivity), which in the absence of refractoriness would only support travelling fronts. We construct these solutions numerically as stationary periodic solutions in a co-moving frame (of both an equivalent delay differential model as well as the original delay integro-differential model). Continuation methods are used to obtain the dispersion curve for periodic travelling waves (speed as a function of period), and found to be reminiscent of those for spatially extended models of excitable tissue. A kinematic analysis (based on the dispersion curve) predicts the onset of wave instabilities, which are confirmed numerically.
Explicit computational formulas for coefficients of the periodic normal forms of the three most complex codim 2 bifurcations of limit cycles with dimension of the center manifold equal to 4 or to 5 in generic autonomous ODEs are derived. The resulting formulas are independent of the dimension of the phase space and involve solutions of certain boundary-value problems as well as multilinear functions from the Taylor expansion of the ODE right-hand side near the cycle. The formulas allow one to distinguish between the complicated bifurcation scenarios which can happen near these codim 2 bifurcations of limit cycles, where 3-tori and 4-tori can be present. We apply our techniques to the study of a known laser model, a novel model from population biology, and one for mechanical vibrations; these models exhibit Limit Point--Neimark-Sacker, Period-Doubling--Neimark-Sacker and double Neimark-Sacker bifurcations. Lyapunov exponents are computed to numerically confirm the results of the normal form analysis, in particular with respect to the existence of stable invariant tori of various dimensions and chaos.
The EEG of patients in non-convulsive status epilepticus (NCSE) often displays delta oscillations or generalized spike-wave discharges. In some patients, these delta oscillations coexist with intermittent epileptic spikes. In this study we verify the prediction of a computational model of the thalamo-cortical system that these spikes are phase-locked to the delta oscillations. We subsequently describe the physiological mechanism underlying this observation as suggested by the model. It is suggested that the spikes reflect inhibitory stochastic fluctuations in the input to thalamo-cortical relay neurons and phase-locking is a consequence of differential excitability of relay neurons over the delta cycle. Further analysis shows that the observed phase-locking can be regarded as a stochastic precursor of generalized spike-wave discharges. This study thus provides an explanation of intermittent spikes during delta oscillations in NCSE and might be generalized to other encephathologies in which delta activity can be observed.
Objective: Characterization of the functional neuronal activity and connectivity within the subthalamic nucleus
(STN) in patients with Parkinson's disease (PD).
Methods: Single units were extracted from micro-electrode recording (MER) of 18 PD patients who underwent STN deep brain stimulation (DBS) surgery. The firing rate and pattern of simultaneously recorded spike trains and their coherence were analyzed. To provide a precise functional assignment of position to the observed activities, for each patient we mapped its classified multichannel STN MERs to a generic atlas representation with a sensorimotor part and a remaining part.
Results: Within the sensorimotor part we found significantly higher mean firing rate (P < 0.05) and significantly more burst-like activity (P < 0.05) than within the remaining part. The proportion of significant coherence in the beta band (13-30 Hz) is significantly higher in the sensorimotor part of the STN than elsewhere (P=0.015).
Conclusions: The STN sensorimotor part distinguishes itself from the remaining part with respect to beta coherence, firing rate and burst-like activity and postoperatively was found as the preferred target area. Significance: Our firing behavior analysis may help to discriminate the STN sensorimotor part for the placement of the DBS electrode.
A lumped model of neural activity in neocortex is studied to identify regions of multi-stability of both steady states and periodic solutions. Presence of both steady states and periodic solutions is considered to correspond with epileptogenesis. The model, which consists of two delay differential equations with two fixed time lags, is mainly studied for its dependency on varying connection strength between populations. Equilibria are identified, and using linear stability analysis, all transitions are determined under which both trivial and non-trivial fixed points lose stability. Periodic solutions arising at some of these bifurcations are numerically studied with a two-parameter bifurcation analysis.
The purpose of this article to provide an explicit example where continuation based on the homotopy method fails. The example is a one-parameter homotopy for periodic orbits between two well-known nonlinear systems, the normal form of the Hopf bifurcation and the van der Pol system. Our analysis shows that various types of obstructions can make approximation over the whole range of the homotopy parameter impossible. The Hopf-van der Pol system demonstrates that homotopy methods may fail even for seemingly innocent systems.
We present a computational model of a thalamocortical relay neuron for exploring basal ganglia thalamocortical loop behavior in relation to Parkinson's disease and deep brain stimulation (DBS). Previous microelectrode, single-unit recording studies demonstrated that oscillatory interaction within and between basal ganglia nuclei is very often accompanied by synchronization at parkinsonian rest tremor frequencies (3-10 Hz). These oscillations have a profound influence on thalamic projections and impair the thalamic relaying of cortical input by generating rebound action potentials. Our model describes convergent inhibitory input received from basal ganglia by the thalamocortical cells based on characteristics of normal activity, and/or low-frequency oscillations (activity associated with Parkinson's disease). In addition to simulated input, we also used microelectrode recordings as inputs for the model. In the resting state, and without additional sensorimotor input, pathological rebound activity is generated for even mild parkinsonian input. We have found a specific stimulation window of amplitudes and frequencies for periodic input, which corresponds to highfrequency DBS, and which also suppresses rebound activity for mild and even more prominent parkinsonian input. When low-frequency pathological rebound activity disables the thalamocortical cell's ability to relay excitatory cortical input, a stimulation signal with parameter settings corresponding to our stimulation window can restore the thalamocortical cell's relay functionality.
The pedunculopontine nucleus has been suggested as a target for DBS. In this paper we propose a single compartment computational model for a PPN Type I cell and compare its dynamic behavior with experimental data. The model shows bursts after a period of hyperpolarization and spontaneous firing at 8 Hz. Bifurcation analysis of the single PPN cell shows bistability of fast and slow spiking solutions for a range of applied currents. A network model for STN, GPe and GPi produces basal ganglia output that is used as input for the PPN cell. The conductances for projections from the STN and the GPi to the PPN are determined from experimental data. The resulting behavior of the PPN cell is studied under normal and Parkinsonian conditions of the basal ganglia network. The effect of high frequency stimulation of the STN is considered as well as the effect of combined high frequency stimulation of the STN and the PPN at various frequencies. The relay properties of the PPN cell demonstrate that the combined high frequency stimulation of STN and low frequency (10 Hz, 25 Hz, 40 Hz) stimulation of PPN hardly improves the effect of exclusive STN stimulation. Moreover, PPN-DBS at low stimulation amplitude has a better effect than at higher stimulation amplitude. The effect of PPN output on the basal ganglia is investigated, in particular the effect of STN-DBS and/or PPN-DBS on the pathological firing pattern of STN and GPe cells. PPN-DBS eliminates the pathological firing pattern of STN and GPe cells, whereas STN-DBS and combined STN-DBS and PPN-DBS eliminate the pathological firing pattern only from STN cells.
Two models of the neocortex are developed to study normal and pathologic neuronal activity. One model contains a detailed description of a neocortical microcolumn represented by 656 neurons, including superficial and deep pyramidal cells, four types of inhibitory neurons, and realistic synaptic contacts. Simulations show that neurons of a given type exhibit similar, synchronized behavior in this detailed model. This observation is captured by a population model that describes the activity of large neuronal populations with two differential equations with two delays. Both models appear to have similar sensitivity to variations of total network excitation. Analysis of the population model reveals the presence of multistability, which was also observed in various simulations of the detailed model.
In this computational study, we investigated (i) the functional importance of correlated basal ganglia (BG) activity associated with Parkinson's disease (PD) motor symptoms by analysing the effects of globus pallidus internum (GPi) bursting frequency and synchrony on a thalamocortical (TC) relay neuron, which received GABAergic projections from this nucleus; (ii) the effects of subthalamic nucleus (STN) deep brain stimulation (DBS) on the response of the TC relay neuron to synchronized GPi oscillations; and (iii) the functional basis of the inverse relationship that has been reported between DBS frequency and stimulus amplitude, required to alleviate PD motor symptoms [A. L. Benabid et al. (1991)Lancet, 337, 403-406]. The TC relay neuron selectively responded to and relayed synchronized GPi inputs bursting at a frequency located in the range 2-25 Hz. Input selectivity of the TC relay neuron is dictated by low-threshold calcium current dynamics and passive membrane properties of the neuron. STN-DBS prevented the TC relay neuron from relaying synchronized GPi oscillations to cortex. Our model indicates that DBS alters BG output and input selectivity of the TC relay neuron, providing an explanation for the clinically observed inverse relationship between DBS frequency and stimulus amplitude.
Persistence and bifurcations of Lyapunov manifolds can be studied by a combination of averaging-normalization and numerical bifurcation methods. This can be extended to infinite-dimensional cases when using suitable averaging theorems. The theory is applied to the case of a parametrically excited wave equation. We find fast dynamics in a finite, resonant part of the spectrum and slow dynamics elsewhere. The resonant part corresponds with an almost-invariant manifold and displays bifurcations into a wide variety of phenomena among which are 2- and 3-tori.
The theory of dynamical systems studies the behavior of solutions of systems, like nonlinear ordinary differential equations (ODEs), depending upon parameters. Using qualitative methods of bifurcation theory, the behavior of the system is characterized for various parameter combinations. In particular, the catalog of system behaviors showing qualitative differences can be identified, together with the regions in parameter space where the different behaviors occur. Bifurcations delimit such regions. Symbolic and analytical approaches are in general infeasible, but numerical bifurcation analysis is a powerful tool that aids in the understanding of a nonlinear system. When computing power became widely available, algorithms for this type of analysis matured and the first codes were developed. With the development of suitable algorithms, the advancement in the qualitative theory has found its way into several software projects evolving over time. The availability of software packages allows scientists to study and adjust their models and to draw conclusions about their dynamics.
We present new or improved methods to continue heteroclinic and homoclinic orbits to fixed points in iterated maps and to compute their fold bifurcation curves, corresponding to the tangency of the invariant manifolds. The proposed methods are applicable to general n-dimensional maps and are implemented in matlab. They are based on the continuation of invariant subspaces (CIS) algorithm, which is presented in a novel way. The systems of defining equations include the Riccati equations appearing in CIS for bases of the generalized stable and unstable eigenspaces. We use the bordering techniques to continue the folds, and provide full algorithmic details on how to treat the Jacobian matrix of the defining system as a sparse matrix in matlab. For a special -- but important in applications -- case n=2 we describe the first matlab implementation of known algorithms to grow one-dimensional stable and unstable manifolds of the fixed points of noninvertible maps. The methods are applied to study heteroclinic and homoclinic connections in the generalized Hénon map.
The paper provides full algorithmic details on switching to the continuation of all possible codim 1 cycle bifurcations from generic codim 2 equilibrium bifurcation points in n-dimensional ODEs. We discuss the implementation and the performance of the algorithm in several examples, including an extended Lorenz-84 model and a laser system.
Bifurcation software is an essential tool in the study of dynamical systems. From the beginning (the first packages were written in the 1970's) it was also used in the modelling process, in particular to determine the values of critical parameters. More recently, it is used in a systematic way in the design of dynamical models and to determine which parameters are relevant. MatCont and Cl_MatCont are freely available matlab numerical continuation packages for the interactive study of dynamical systems and bifurcations. MatCont is the GUI-version, Cl_MatCont is the command-line version. The work started in 2000 and the first publications appeared in 2003. Since that time many new functionalities were added. Some of these are fairly simple but were never before implemented in continuation codes, e.g. Poincaré maps. Others were only available as toolboxes that can be used by experts, e.g. continuation of homoclinic orbits. Several others were never implemented at all, such as periodic normal forms for codimension 1 bifurcations of limit cycles, normal forms for codimension 2 bifurcations of equilibria, detection of codimension 2 bifurcations of limit cycles, automatic computation of phase response curves and their derivatives, continuation of branch points of equilibria and limit cycles. New numerical algorithms for these computations have been published or will appear elsewhere; here we restrict to their software implementation. We further discuss software issues that are in practice important for many users, e.g. how to define a new system starting from an existing one, how to import and export data, system descriptions, and computed results.
We discuss new and improved algorithms for the bifurcation analysis of fixed points and periodic orbits (cycles) of maps and their implementation in matcont, a MATLAB toolbox for continuation and bifurcation analysis of dynamical systems. This includes the numerical continuation of fixed points of iterates of the map with one control parameter, detecting and locating their bifurcation points (i.e., limit point, period-doubling, and Neimark-Sacker) and their continuation in two control parameters, as well as detection and location of all codimension 2 bifurcation points on the corresponding curves. For all bifurcations of codim 1 and 2, the critical normal form coefficients are computed, both numerically with finite directional differences and using symbolic derivatives of the original map. Using a parameter-dependent center manifold reduction, explicit asymptotics are derived for bifurcation curves of double and quadruple period cycles rooted at codim 2 points of cycles with arbitrary period. These asymptotics are implemented into the software and allow one to switch at codim 2 points to the continuation of the double and quadruple period bifurcations. We provide two examples illustrating the developed techniques: a generalized Henon map and a juvenile/adult competition model from mathematical biology.
We study codimension-2 bifurcations of fixed points of dissipative diffeomorphisms with a pair of complex eigenvalues together with either an eigenvalue $-1$ or another such a pair. In the previous studies only cubic normal forms were considered. However, in some cases the unfolding requires higher order terms and these are investigated here. We (re)derive the normal forms and reduce them to a single amplitude map. This map is similar to the amplitude system for the double-Hopf bifurcation for vector fields. We show how the critical normal form coefficients determine the general bifurcation picture for this amplitude map. Representative nonsymmetric perturbations of the normal forms are studied numerically. Our case studies show a detailed picture near various bifurcation curves, which is richer than known theoretical predictions. For arbitrary maps with these bifurcations we give explicit formulas for critical normal form coefficients on center manifolds. We apply them to an example from robotics where we are able to demonstrate the existence of a bubble-structure, which was only observed in perturbed normal forms before.
The Generalized Hénon Map(GHM) appears in codim 2 bifurcation of homoclinic orbits for maps, for example a neutral saddle with a homoclinic tangency (NHT). The standard Hénon map with the standard parameters itself is interesting for the strange attractor, the homoclinic orbit and the resulting chaos. However, from the bifurcation point of view it turns out that many appearing bifurcations of the standard Hénon map are degenerate, for example, the strong resonances. This is not the case for GHM. One interesting feature is that the GHM has the NHT-bifurcation, so we show that the existence of infinitely many of such bifurcations.
We derive the formulas for the center manifold reduction and the critical normal form coefficients for the codimension two bifurcations for diffeomorphisms of the plane. Some formulas involve third or fifth order derivatives. This poses no problems for maps, but for limit cycles of ODE's it does. A RungeKuttaFehlberg-78 integrator adapted with automatic differentiation(ADOL-C 1.8) is used to obtain the required data.
The fold-flip bifurcation occurs if a map has a fixed point with multipliers +1 and -1 simultaneously. In this paper the normal form of this singularity is calculated explicitly. Both local and global bifurcations of the unfolding are analysed by exploring a close relationship between the derived normal form and the truncated amplitude system for the fold-Hopf bifurcation of ODEs. Two examples are presented, the generalized H\'enon map and an extension of the Lorenz-84 model. In the latter example the first, second and third order derivatives of the Poincar\'e map are computed to find the normal form coefficients.
We study a model for a robot-arm due to Steindl, Troger and Lindtner near a bifurcation of a fixed point which has multipliers +1 and -1 simultaneously. There is a degeneracy in the parameter dependent part, so that the multiplier +1 is associated to a transcritical bifurcation instead of an ordinary fold. The critical center-manifold reduction is given and the nonstandard unfolding as well.