Introduction to Focus Issue: Dynamical disease: A translational approach

Intrinsic to biological systems in general, and physiological systems in particular, are fluctuations in time of their state variables. These variations may take place over long time intervals and may be hard to identify, their clinical relevance notwithstanding. Recognizing them in his own clinical practice, Reimann² documented more than 2000 clinical cases of what he called “periodic diseases,” various medical disorders recurring at time intervals longer than 24 h. Many of these cases concerned hormonal³ and hematological^4,5 regulations. A seminal breakthrough occurred when chaotic behavior was observed in a delay differential equation (DDE) modeling respiratory and hematopoietic systems⁶ in what is now called the Mackey–Glass equation (but could have been known, had history been slightly perturbed,⁷ as the Lasota equation⁸).

This equation was instrumental in testing a number of techniques for dynamical time series analysis, following an innovative phase-plane construction for infinite-dimensional systems.⁹ In the following decade, nonlinear dynamics was explored for its relevance in physiological systems, in both normal and pathological conditions.¹⁰ This eventually led to the workshop “Dynamical Disease” in February 1994, which specifically addressed the use of nonlinear dynamics technique in a clinical setting and the ensuing Focus Issue of CHAOS.¹ Many of the physiological systems investigated then, and their pathologies, are still the object of considerable interest. The current state of the art in the understanding of many of these dysregulations in terms of bifurcations, and the extent to which it is the relevant approach, are presented in the current Focus Issue.

The challenge of translating mathematical advances into medical procedures is a never ending endeavor.¹¹ The traditional approach of proposing mathematical models for specific physiological systems and carrying out analysis of the properties of these models using nonlinear dynamics has been implemented extensively. However, the next step of translating mathematical insights into practical medical applications at the bedside remains far from straightforward. The recent explosion of data-collecting devices, software, and environments is rapidly changing the way that health care professionals diagnose and monitor diseases. Thus, the time seems ripe to try to reconcile these “hypothesis-driven” and “data-driven” approaches.

This Focus Issue comprises Refs. 14, 19, 20, 24–26, 28–31, 34–36, 38, 39, 40, 44, 47, 48, 55, 56, 59, 62, 73, and 76. Historical perspective on the influence of the seminal paper⁶ and its authors has also recently appeared.^12,13

Mackey¹⁴ evaluated the concept of a dynamical disease from the perspective of periodic hematological diseases. Although these diseases are rare, they have been modeled extensively as examples of dynamical diseases (for a recent review, see Ref. 15). The quintessential example of a dynamical disease is cyclic neutropenia (CN). The control mechanism appears to be functioning normally; however, it operates in an “abnormal” parameter range in which there is increased apoptosis (i.e., programmed cell death) in the proliferating neutrophil precursors, but not in the neutrophil stem cells. It should be noted that the mathematical models for CN predict the possibility that the locally stable limit cycle may co-exist with a locally stable steady state. This has been observed for gray collies, an animal model for CN,¹⁶ but not yet confirmed for the human disorder. The dynamical disease status of the other periodic blood cell disorders remains very much a work in progress: periodic chronic myelogenous leukemia (maybe yes), cyclic thrombocytopenia (unclear), and periodic autoimmune hemolytic anemia (unclear).

A difficulty with the original Mackey–Glass definition of a dynamical disease is understanding what is meant by the expression “a basically normal control system.” Surely it is likely true that patients in whom control systems have been altered by disease would also experience changes in dynamical behaviors that arise because of parameter changes. Moreover, changes in dynamics do not always reflect changes in parameters but can arise, for example, because of the effects of perturbations and noisy stimuli. In other words, a strict application of the definition of a dynamical disease would eliminate most of the articles that appear in this special issue on dynamical diseases! It has been suggested that the term “dynamic disease” would be more inclusive.^17,18

A case in point concerns the onset patterns of the focal-onset epilepsies. The location of seizure onset, i.e., the epileptic focus, is always associated with anatomical and functional abnormalities. A wide range of complex patterns are observed in the electroencephalogram (EEG) at seizure onset including sharp activity, spike-and-wave discharges, bursts of high amplitude polyspikes, burst suppression, and delta brushes. However, the relationship between the EEG patterns and the underlying anatomical and functional abnormalities that underlie these seizure onset patterns has not yet been elucidated. Zhang et al.¹⁹ argue that the different onset patterns can be explained by a slow shift of a bifurcation parameter in and out of a range where an oscillatory state (seizure) occurs. They first demonstrate this mechanism in a three-variable model in which an oscillator is described by a two variable model that is connected to a very slow population in a feedback loop. Provided that the parameters in the equations that describe the oscillator are tuned near a suitable bifurcation point, spontaneous transitions between a non-oscillatory and oscillatory state occur. This model is then extended to a minimal neural network in which two connected cortical oscillators are under the control of one common slow time scale population. All of the observed clinical patterns for focal seizure onset could be generated by fine tuning the parameters in the model. Although the nature of the slowly changing parameter is not known, one possibility is that it might be related to changes in the extracellular space regulated by neural–glial cell interactions.

A second example arises in the setting of models for cancer. Balti et al.²⁰ illustrate the use of dynamical systems analysis to provide insights into the behaviors of a well established model that describes the interaction between tumor and immune system in response to the administration of radiotherapy and immunotherapy.²¹ This model takes the form of a six-dimensional system of nonlinear ordinary differential equations (ODEs). The T cell tumor tissue infiltration rate, $SL$⁠, is the important factor responsible for individual therapeutic inter-variability. The authors show that changing $SL$ can produce a saddle-node bifurcation point which causes this dynamical system to oscillate between two states that correspond, respectively, to the tumor-free and maximum tumor volume states. When $SL$ is low, these two states co-exist, whereas when $SL$ is high, only the maximum tumor volume state exists. The performance of therapeutic strategies is discussed from the point of view of the nature of the basins of attraction of these two states. In particular, if therapy is considered to be the input, then this translates into changes in the trajectory shapes of the solutions as they approach the equilibrium points.

Dendritic cells (DCs) play a critical role in linking adaptive and innate immune responses. Although it seems likely that dendritic cells would be important in immunotherapeutic protocols for the treatment of cancer, clinical trials so far have been disappointing. Current mathematical models that describe the dynamics of tumor-immune interactions in the presence of dendritic cells take the form of a four-dimensional system of ordinary differential equations (ODEs).^22,23 Dickman and Kuang²⁴ examine the effect of including a time delay into these models. The time delay accounts for the time it takes for activated cytotoxic T lymphocytes to become effective in killing tumor cells. Using clinical observations, the authors are able to reduce the ODE models to a second-order delay differential equation (DDE), by first assuming a constant proportion of DCs in the blood and then applying a quasi-steady state approximation for them, a simplification justified by clinical observations. This second-order DDE exhibits complex dynamics including a Hopf bifurcation, a reverse Hopf bifurcation, and bistability. For small delays, the tumor cells and immune cells co-exist, whereas large delays produce fatally high tumor cell levels. Interestingly, for intermediate values of the time delay, very complex periodic and even chaotic oscillations can occur. In this regime, the outcomes are uncertain and depend on the choice of initial tumor cell burdens. It is not difficult to imagine that the chaotic nature of the responses in this regime would confound the outcome of clinical trials to test the effectiveness of DCs. The model suggests a number of strategies that might be used to mitigate these effects.

What is the mechanism that makes it possible for cancer cells to invade healthy tissues? One possibility is that it is related to the dynamics exhibited by “active matter” in confinement. The term active matter refers to molecules that consume energy to generate directed motion. Zhao et al.²⁵ investigate the stability of the shear flow of active extensile filaments which are confined within a narrow channel. The self-propelled elongated fiber units take the form of bundles of microtubules powered by an adenosine triphosphate (ATP)-consuming kinesin 23 molecule attached at one end. The mathematical model takes the form of the Erickson–Leslie equations of the balance of linear and angular momentum. The role of kinesin 23 is modeled by adding an activity source term to the Leslie stress. It is shown that the stability of the flow depends upon the geometry and magnitude of the flow profiles, a threshold value of the activity parameter, and the channel aspect ratio. The generalized eigenvalue problems that arise in the stability analysis are solved using a spectral method, which incorporates Chebyshev polynomials.

A formidable challenge for obtaining mathematical insights into clinical observations is the need to estimate the underlying model parameters. Often, investigators “guess” reasonable values of the parameters (“guestimate”) and then perform a sensitivity analysis and/or draw the parameter values from a statistical distribution in order to estimate the effects of a wrong estimate. However, it is much better to estimate the parameters experimentally. Two new methods for parameter estimation were discussed.

Postural sway refers to the fluctuation in the center of pressure measured while a subject stands quietly on a force platform with eyes closed. Measurements of postural sway are important for identifying elderly subjects who are most at risk of falling. Surprisingly, the controlling forces for maintaining balance are exerted intermittently.^26–28 Suzuki et al.²⁸ introduce the technique of parameter inference using approximate Bayesian computation for the estimation of parameters in a mathematical model for intermittent control of balance. This method provides a framework for identifying parameter values in the model, which produce data that are sufficiently similar to the observed data as evaluated by a likelihood function. Their intermittent control model is parameterized by a parameter that is associated with the degree of intermittency. Thus, it is possible to identify individuals who use an intermittent balance control strategy and those that do not. They apply this technique to healthy individuals and those with Parkinson’s disease (PD). It is observed that healthy individuals and subjects with mild postural symptoms use an intermittent feedback control strategy for maintaining balance. In contrast, the elderly and patients with severe postural symptoms use a continuous feedback control strategy.

Bahari et al.²⁹ introduce an iterative prediction–correction scheme referred to as unscented Kalman filter based data assimilation (DA) to estimate measured and “hidden” variables in a manner that is capable of predicting the future state of the system. This method iteratively couples the mathematical equations of a computational model to ongoing experimental measurements of the system. It is shown that once the error between the tracked and its corresponding true state is minimized, it is possible to reliably reconstruct the unobserved system variables. Thus, the initial computational model may not need to be sophisticated. This procedure is illustrated through a study of sleep-wake dynamics in freely behaving rodents. Even though the starting mathematical model does not possess a circadian drive, DA correctly estimated the 24-h sleep-wake cycles as well as the underlying interaction with the suprachiasmatic nucleus!

The relationships between dose, exposure, and drug effects complicate the understanding of the pathophysiological aspects traditionally addressed in mathematical pharmacology. Are the abnormal dynamics observed at the bedside the result of the pathophysiology of the disease and/or the result of the physicians attempts to treat the disease with pharmacological agents? In other words, mathematical models for the treatment of a disease need to consider both what the drug does to the body (pharmacodynamics) and what the body does to the drug (pharmacokinetics). These considerations have motivated the development of quantitative systems pharmacology modeling techniques to understand the disease–treatment interactions.

Ursino et al.³⁰ and Véronneau-Veilleux et al.³¹ team up to provide a provocative illustration of this approach by examining how the effects of levodopa for the treatment of Parkinson’s disease (PD) change as the disease severity increases. The model has two parts: a system of nonlinear ODEs that describe a two-compartment model for the pharmacokinetics of levodopa and a neurocomputational model that describes activity in the basal ganglia (BG),^32,33 i.e., the site of action of levodopa. Nonlinearities in dopamine dynamics arise because of Michaelis–Menton enzyme kinetics whereas the nonlinearities in neural networks arise from the sigmoidal activation functions of neurons. It is shown that the pharmacokinetics of levodopa depend on the degree of nigrostriatal denervation, $f$⁠. A consequence is that a plot of dopamine concentration in the plasma and brain exhibits hysteresis; in other words, there is a delay between when the drug peaks in the plasma and when it peaks in the brain. The effects of this hysteresis on the patient’s motor responses are examined using a nonlinear deterministic neurocomputational model of the BG that describes finger tapping, a clinically relevant output. This model describes the role of the BG to select actions and hence includes the direct or Go pathways, the indirect or NoGo pathways, and the hyperdirect pathways via the action of the sub-thalamic nucleus. Depending on the model pathways, this model can describe a wide variety of movement abnormalities seen in patients with PD. To model finger tapping, the Go pathways are identified with finger up and the NoGo pathways with finger down. During mild Parkinson’s disease (⁠$f∈[0.1,1]$⁠), finger tapping can remain high for a long time even when the levodopa plasma concentration is low. However, for severe disease (⁠$f<0.1$⁠), finger tapping first increases and then exhibits a dramatic fall. These complexities resemble the complexities encountered by clinicians and patients as they manage their symptoms.

Health care professionals are often required to make difficult decisions based on limited data. For example, what is the expected survival time of their patients (i.e., risk stratification, how can the success of a treatment be assessed, and so on)? Can mathematical models help guide these professionals in making these decisions?

An important step in the clinical management of patients with cancer is determining the stage of a patient’s disease. In other words, is the cancer localized or has it spread to the rest of the body, does the patient test positive to certain biomarkers known to affect survival and response to therapy, and so on. Typically, risks are assigned empirically from the clinical course, radiological findings, and the presence of a variety of molecular and cytogenic markers. Determining the patient’s prognosis makes it possible to adapt therapy to the individual needs of the patient. Stiehl³⁴ shows how a mathematical model of stem cell dynamics can be used to validate the use of counts of blood-producing stem cells (HSCs) as a biomarker for risk stratification in patients who suffer from acute myeloid leukemia. The basic idea is that leukemic stem cells (LSCs) out-compete the healthy HSC in a stem cell niche that has fixed finite capacity. Thus, lower HSCs correlate with more severe disease. The mathematical model takes the form of 12 nonlinear ordinary differential equations. It is shown that the model has three non-negative steady states: a trivial steady state, a healthy steady state, and a fully leukemic steady state. A numerical sensitivity analysis indicates that the LC proliferation rate and HSC dislodgement probability are the key determinants of the speed of disease progression and that the parameters that describe the leukemic non-stem cells have a negligible effect. The dynamics strongly support the empirical suggestion that HSC counts are a suitable surrogate marker for both LSC counts and kinetics.

Surgical removal of the epileptic focus is a potentially life-changing treatment for patients whose seizures cannot be adequately controlled by anticonvulsant medications. Over the last 20 years, it has become established that focal-onset seizures can be produced by a widespread cortical network.¹⁸ This may explain why in some patients seizures can re-appear after surgery. Junges et al.³⁵ examine the consequences of node removal, a surrogate for surgery, on the dynamics of the remaining network of excitable nodes. Each node is modeled by a modified version of the normal form for a subcritical Hopf bifurcation. An additional equation is included to describe the slow variations of an excitability variable. There are conduction time delays between the nodes. The authors consider the 199 non-isomorphic networks composed of four nodes that are weakly or strongly connected and use numerical simulations to evaluate the effects of node removal. They define the Brain Network Ictogenicity (BNI) as the proportion of time that nodes within the network spend in a seizure-like state. In general, the BNI decreases as the number of edges increases. However, the effects of the number of edges also depend on the topology of the network and the hierarchy of the connections. A consequence is that networks with a high BNI could be identified in which the BNI could be reduced by the selective removal of a specific node. However, in these cases, any alteration to the remaining network structure resulted in an increase in the BNI, even to values as high as of those observed prior to node removal! A similar observation was obtained for approximately $10%$ of networks composed of ten nodes. These observations suggest the importance of consideration of the dynamics of epileptic networks in the design of surgical resections in the treatment of epilepsy.

The mechanisms involved in human disease and its treatment are exceedingly complex. This complexity is a consequence of the nonlinear nature of the underlying dynamical systems. Proposed therapies must run through a gauntlet of clinical trials before they are approved. A fundamental problem is that potential pitfalls as well as the identification of optimal treatment regimens are not known a priori. The consequent inefficiency of this process translates into the high cost of medical care. Alfonso et al.³⁶ review the current status of in silico clinical trials for the treatment of certain cancers and systemic inflammation. An in silico clinical trial involves a mathematical model that describes the disease process as well as the pharmacokinetics and pharmacodynamics of the drug. Virtual patients are created by randomly drawing the values of important parameters from a plausible distribution. The drug trial is performed using computer simulations, hence the term in silico. The results to date emphasize the potential for in silico drug trials to optimize treatment protocols and identify pitfalls. However, the article does not indicate the current status of these approaches in the design of ongoing drug trials. It can be anticipated that the use of in silico drug trials will have medico-legal implications and also impact curricula for medical education.

The nature of health care is dramatically changing. The driving force is the powerful computational capabilities of the personal cell phone. It has been estimated that Americans spend 3–5 h per day using their cell phone and this usage time is ever increasing. Thus, acting together with developments in sensors, cell phone applications can not only continuously monitor key physiological variables (temperature, activity, blood sugar), but automatically alert healthcare workers when attention is necessary. It is not difficult to imagine that soon cell phone applications will become an essential component of health management. For example, continuous monitoring of movement in patients with PD³⁷ will make it possible to tailor pharmacotherapies to individual patient needs and provide a patient-centrist approach to individualized medical care.

Cell phone applications that record the electrocardiogram (ECG) of the heart for long periods of time (days, weeks, years) are already available. Glass³⁸ outlines two problems related to heart disease that could be aided by cell phone ECG monitoring: (1) diagnosis of atrial fibrillation based on the time intervals between subsequent beats, and (2) risk stratification for cardiac sudden death based on the nonlinear dynamic mechanisms of premature ventricular complexes (PVCs) (see also Ref. 39). Indeed, implantable cardiac defibrillators that monitor the cardiac rhythm and then deliver a shock to abort a life-threatening cardiac arrhythmia are already regarded as the standard of care for patients at risk. The problem is to identify those patients who require implantation. It is possible that understanding the long-term dynamics of PVCs may help identify patients who require implantation.

Twose et al.⁴⁰ provide an important “proof of concept” that demonstrates the potential for cell phone applications to identify patients who are having changes in the state of the disease and hence require medical follow-up. Their goal was to use cell phones to detect early warning signs (EWS) of disease worsening in a group of patients with multiple sclerosis (MS). Their assumption was that different people have different typing mannerisms that can be used to make a “typing signature” that is unique for each user. A Neurokeys App is used to collect data such as the differences in the times between key presses and release, the time to correct a typing error, etc. The data streams are found to exhibit many of the properties of a complex dynamical system such as memory, non-stationarity, and sensitive dependence on initial conditions. Critical transitions in the keystroke dynamics are observed for patients with MS but not in healthy controls. Examples of EWS detected by changes in keystroke dynamics are given, which coincide with clinical relevant changes in outcome measures that are found in patients with MS with and without differences in the amount of enhancing lesions in brain magnetic resonance imaging (MRI). This is not unexpected since a worsening of symptoms in patients with multiple sclerosis is not necessarily a sign of active disease. It can be due, for example, to an increase in the ambient room temperature.

The increasing availability of long time series of physiological variables, such as heart rate, cell numbers, and electrical discharges, raises a number of practical problems that require the attention of mathematicians and statisticians. These issues include, but are not limited to, the identification of the important physiologically accessible parameters to monitor, the development of efficient data mining techniques to detect an abnormality, and statistical analytic techniques to identify artifacts and determine levels of significance that would motivate medical intervention.

Physiological time series typically exhibit non-stationarity. The standard approach for the analysis of such time series is to use techniques to remove trends so that the resulting time series becomes stationary! However, when the time series is generated by a nonlinear dynamical system, this approach runs the risk of obscuring information that may provide important insights. A case in point is provided by the occurrence of premature ventricular contractions (PVCs) recorded by the electrocardiogram (ECG). Bury et al.³⁹ examine the patterns of PVCs occurrences in continuous ECG recordings of at least 6 days in length for patients with frequent ectopy. Since the number of known mechanisms that generate PVCs is small, the hope was that the observed dynamics could be classified in terms of the bifurcations generated by mathematical models. For each patient, the pattern of PVCs could be constant over a period of a single day, and well described by an appropriate mathematical model. However, the patterns of PVCs were typically varying from one day to the next. For some patients, the pattern of PVCs depended on their heart rate (HR); for others, it depended on both HR and time of day. It is not difficult to imagine that the slow variations in the PVC patterns might, in some sense, map out the state space for the mechanisms that generate PVCs in a given patient. Thus, this rich dynamical behavior would be expected to provide important patient-specific insights that could perhaps identify those most at risk for a heart attack.

A pressing question for both patients and physicians is whether or not adverse medical events, such as the occurrence of an epileptic seizure, can be forecasted. For example, patients with epilepsy often keep a diary of their seizure occurrences in an effort to search for a pattern.⁴¹ Short-term (⁠$<2$ weeks⁴²) and long-term (up to 10 years⁴³) recordings with implanted electo-encephalogram (EEG) electrodes indicate a degree of predictability and, in particular, suggest a link between seizure occurrence and endogenous circadian (⁠$≈24$ h), multidien (weeks to months) and even circannual (years) cycles. Leguia et al.⁴⁴ ask the question whether it is better to forecast seizure occurrence by fitting cycles of fixed-length to the time series (M1) or by deriving continuous cycles from measurements of a biomarker such as epileptic spikes (M2). Not surprisingly, their results demonstrate that only time frequency analysis of continuous recordings of a related biomarker can reveal the full dynamics of the manifested cyclical behaviors. In particular, the M2 approach provides stronger support for the evidence of a cycle and in addition is more robust in the face of non-stationarities. One explanation of the superiority of the M2 over the M1 methods is that biological rhythms frequently can decouple from external influences. The long-term recordings in this study were made possible because patients were implanted with the Neuropace RNS© system, which both detects the occurrence of a seizure and then delivers a stimulus to abort it.¹⁸ The effect of this treatment strategy on the pattern of seizure recurrences has not been investigated. However, it is quite likely that co-adaptation to implanted brain-machine interfaces (BMIs) occurs.⁴⁵ Dynamic game theoretic approaches raise the possibility that it may be possible to exploit these interactions to optimize the performance of the electroceutical device.⁴⁶ 

The observation that biological oscillators can often be stopped by applying a stimulus of the right size, at the right time, has fueled the development of electroceutical devices for the treatment of, for example, certain cardiac arrhythmias and epileptic seizures.

It is often difficult at the bedside to determine the timing of stimulus to effect a change. Chang et al.^47,48 show that it is possible to construct a stimulus that stops an oscillation that is “phase agnostic,” namely, it is successful no matter at what phase it is applied. The agnostic stimulus, typically a complex waveform, can be determined using an extrema distortion algorithm.⁴⁹ Numerical simulations show that the application of the phase agnostic stimulus first confines, or “corrals,” the oscillation to a narrow range of new phases and then forces the oscillation to its phase singularity. A similar “phase corralling-drive to singularity” scenario is observed for a variety of oscillators including the radial isochron clock, the FitzHugh–Nagumo equations, and mathematical models of epileptic seizures,^50,51 and a cardiac model of re-entrant tachycardia⁵² suggesting that it may be generic. The existence of a phase agnostic stimulus for the Glass–Josephson model for a re-entrant tachycardia⁵² is particularly of interest since the oscillations in this model arise both in time and space.

Many diseases of the nervous system are characterized by increased neuronal synchrony. One clinically proven way to decrease neuronal synchrony in patients with Parkinson’s disease is high-frequency deep brain stimulation.⁵³ However, the beneficial effects quickly disappear once the stimulation ceases. A much better way to produce longer lasting desynchronization of neural populations is to stimulate the population at multiple sites using phase-shifted stimuli. The theoretical basis for this approach is called coordinated reset (CR).⁵⁴ Kromer et al.⁵⁵ use computer simulations to determine the optimal number of CR stimulation sites. Specifically, they investigate the effect of CR stimulation on a network of recurrently connected excitatory neurons. The neurons are modeled as integrate-and-fire neurons with spike-timing-dependent-plasticity (STDP). The simulations show that the long-lasting effects of CR stimulation are more pronounced if the stimulation parameters are adjusted with respect to the characteristics of STDP rather than based on the stimulation frequency. There is a nonlinear dependence between the number of stimulation sites and the CR frequency. The fact that the benefits of CR stimulation are not an increasing function of the number of stimulation sites suggests that it might be possible to implement this treatment with a manageable number of stimulation sites.

The organization of large numbers of neurons into synchronized waves lies at the basis of phenomena such as consciousness, cognition, vision, motor behaviors, and the generation of the EEG. It is possible that rehabilitation requires that the normal traveling waves of cortical activity, initiated by localized epicenters of the brain, be re-activated to coordinate neural population activities in both space and time. Budzinskiy et al.⁵⁶ perform a bifurcation analysis of a neural field model proposed by Ref. 57 in order to determine the conditions for traveling wave initiation in an epicenter. This model takes the form of two nonlinear integrodifferential equations with symmetric connectivity functions and no time delay. The analysis demonstrates that depending on the values of the parameters, two types of mutually exclusive solutions can arise: a standing oscillating solution and a periodic traveling wave. The standing wave solutions arise when the local forcing terms in the model are symmetric. The traveling wave solutions arise when the forcing terms are asymmetric. In this case, if the periodic wave is stable, then the solutions converge to the wave. An important observation is that the direction that the traveling wave propagates is also determined by the forcing terms. Indeed in 2D, it is possible to direct traveling waves in any direction between $0$ and $2π$⁠. This observation suggests the potential for epicenters to exchange directed signals. Thus, it may be possible, for example, to use external brain stimulation to restore normal wave propagation to treat chronic situations such as post-stroke aphasia.⁵⁸ 

A naturally occurring situation in which pulsatile control is relevant arises in the setting of neuroendocrine control. Whereas neurons communicate by discrete pulses, the pituitary gland and its target organs interact with continuous feedback. Churilov et al.⁵⁹ introduce a functional–differential equation model for neuroendocrine control, which combines continuous Goodwin-type feedback control with impulsive inputs generated by the hypothalamus using an integrate-and-fire mechanism. In the absence of hypothalamic input, the feedback does not generate oscillations. In the presence of the pulsatile inputs, the model generates a variety of realistic hormonal profiles including ultradian and circadian rhythms and even chaotic dynamics. The results demonstrate that the temporal pattern of the pulsed feedback is no less important than the hormone concentrations is generating the dynamical patterns.

The effects of noisy stimuli on neural dynamics continue to attract attention as a possible therapeutic strategy.^60,61 Powanwe and Longtin⁶² examine the effects of state-dependent (“multiplicative”) and state-independent (“additive”) noise on the generation of rhythms in the beta and gamma range in a stochastic Wilson–Cowan model for cortical dynamics. Acting together, the beta and gamma rhythms generate randomly occurring bursts of higher amplitude oscillations that are thought to be important for neural function. In the absence of noise, the parameters of the Wilson–Cowan model are tuned into a parameter region where the dynamics spiral in to a fixed-point. They use an envelope-phase analysis to show that in the presence of multiplicative noise, the real part of the “effective eigenvalues” becomes less negative as the noise intensity increases. Indeed, for sufficiently intense multiplicative noise, the effective eigenvalue can become positive. A detrimental consequence would be increased synchronization of bursts as seen in certain motor diseases; however, in other contexts, this effect could be beneficial.⁶³ 

Another situation in which the interplay between noise and delay is important is human balance control.²⁸ Falls in aging western societies are a significant cause of morbidity and mortality. Grebrószki et al.²⁶ examine the effects of time delay, sensory quantization, and control torque saturation on stabilizing an inverted pendulum. Both sensory quantization and control torque saturation mitigate the effects of overcontrol, which arises from the interplay between additive noise and time-delayed feedback.⁶⁴ However, quantization does not increase $τcrit$⁠, i.e., the longest time delay, $τ$⁠, for which a pole of length $ℓ$ can be stabilized. Instead, it increases the size of the D-shaped region in the plane of the proportionate and derivative gains for which the upright pendulum position is stable.⁶⁵ The control gains, which are most robust to the effects of multiplicative noise, are located in the lower left hand corner of this D-shaped region. This has also been observed experimentally for pole balancing on the fingertip, virtual stick balancing, and the response of the standing balance to perturbations.

Artificial neural networks (ANNs) are increasingly being used to search for patterns in large amounts of data. Thus, it is not difficult to imagine that ANNs will play an important role in bringing concepts from dynamical diseases to the bedside. Although ANNs are loosely based on the anatomy of the human brain, they do not typically contain time delays. Indeed, the hallmark of the nervous system is the presence of distance-dependent conduction time delays.⁶⁶ However, it is only relatively recently that the study of delayed ANNs has begun to attract attention.^67–71 

It is surprising that the effects of periodic stimulation of a network of neurons are only currently receiving interest. Historically, the ability of the human visual system to detect flickering light was important for the evolving movie industry⁷² and these investigations, in turn, motivated the famous Hartline–Ratcliff equations.⁶⁹ Wang et al.⁷³ investigate the effect of periodic stimulation of a randomly and sparsely connected network of spiking inhibitory and excitatory neurons with time-delayed connectivity.^71,74 Depending on the stimulation frequency, the network responses included resonance, entrainment, oscillation suppression, and a variety of phase locking patterns, some of which have been observed experimentally. Although the authors conclude that stimulation frequency is most important for determining network dynamics, they were able to determine the response function of the network using a mean field approach. In the absence of external stimulation, the response function is nonlinear and resembles a steep sigmoid. Single anodic stimulation and bidirectional phasic stimulation act to linearize the response function similar to the linearization of the response function due to stochastic perturbations.^74,75

Simple neural networks have played an important role in elucidating the roles played by the basal ganglia (BG) in action selection. Cutsuridis et al.⁷⁶ use a neural network to uncover the explanation for the measured difference in the performance of antisaccade tasks between healthy subjects and those with early Huntington’s disease, a neurodegenerative disease that affects the BG. Typically, in response to the sudden appearance of an object in the visual field a subject’s eyes make a saccadic response directed toward the object. In the antisaccade paradigm, the subject is trained not to make a saccade toward the object [error prosaccade (EP) response], but instead to direct the saccade in the opposite direction [correct prosaccade (CP) response]. A one-layer competitive accumulator-to-threshold neural network is constructed: the left half of the network is responsible for generating the EP response and the right half for the CP response. All neurons were connected by short range excitatory and long-range inhibitory connections using a shifted Gaussian kernal, which depended only on the internodal distance. The error decision signal, $Ir$⁠, was inputted by a single node on the left hand side of the network, and the correct decision signal, $Ip$⁠, to a single node in the right hand side of the network. It was assumed that $Ip>Ir$ and that there was a presentation time delay between when the inputs were presented to the network (⁠$Ir$ first). It was shown that poor performance in the antisaccade task can be attributed to the end product of neural lateral inhibition between the CP and EP neural populations. Thus, there is no need to postulate, as assumed previously, that a third “top-down” signal from the cortex interferes with the decision making process.

In a different application, Ursino et al.³⁰ develop a neural network model for action selection in the BG (the direct or GO pathways, the indirect or NoGo pathways and the hyperdirect pathways that act via action of the subthalamic nucleus). Conduction delays are included in the network to account for the delay between one action and the next. This addition was made to enable this model to describe the finger tapping task investigated by Véronneau-Veilleux et al.³¹ (see above). It is shown that this network can produce a variety of dynamical patterns that resemble the movement abnormalities seen in patients with Parkinson’s disease. These patterns include bradykinesia, freezing, repetition, and irregular fluctuations.

The possibility that abrupt changes in state can be predicted has captured the public’s imagination.^77–79 Recent “brain inspired” reservoir computing techniques have taken this possibility one step further by demonstrating that it is possible to accurately forecast complex time series, even those that are chaotic, well into the future.⁸⁰ The ability to predict and forecast makes it possible to make changes in the present that alter future health outcomes. In other words, an impending heart attack or epileptic seizure could be prevented from occurring rather than treating it when it occurs.

It is tempting to imagine that if enough physiological data are collected, methods from artificial intelligence could be used to forecast adverse events. However, it is unlikely that the amount of data required to train such networks will be available due to legal challenges related to the ownership of the data.⁸¹ A surprising limitation of present day approaches that utilize ANNs is that most do not take into account what has already been learned about the dynamical system by the efforts of scientists and mathematicians. We expect that the future will see a merging together of dynamical systems concepts with big data techniques to better understand and treat human illness. Indeed, we may hope to quote from another context, that “we shall see in the next few decades a new generation of mathematical biologists beginning to tackle problems in which the complexity is fundamental.”⁸² 

A number of papers in this Focus Issue were presented at a workshop “Mathematical Physiology: Better Health through Mathematics” held at the Université de Montréal on November 4–22, 2019. The workshop was supported by grants from the Society for Mathematical Biology (SMB) and the Centre de Recherches Mathématiques (CRM), Université de Montréal. J.G.M. acknowledges support from the Simons Foundation and the William R. Kenan Jr. Charitable Trust. F.N. and J.B. acknowledge support from the Natural Sciences and Engineering Research Council (NSERC) of Canada. The authors warmly thank the personnel of the Editorial Office of Chaos, namely, Deborah Doherty, Kristen Overstreet, and Matthew Kershis, for their steady, unrelenting support in making this Focus Issue possible.