# Mackey-Glass equation

Post-publication activity

Curator: Michael Mackey Figure 1: Dynamics in the Mackey-Glass equation, Equation (1), for $$\gamma=1\ ,$$ $$\beta=2\ ,$$ $$\tau=2\ ,$$ and $$n=9.65\ .$$

The Mackey-Glass equation is the nonlinear time delay differential equation $\tag{1} \frac{dx}{dt} = \beta \frac{ x_{\tau} }{1+{x_{\tau}}^n}-\gamma x, \quad \gamma,\beta,n > 0,$

where $$\beta\ ,$$ $$\gamma\ ,$$ $$\tau\ ,$$ $$n$$ are real numbers, and $$x_{\tau}$$ represents the value of the variable $$x$$ at time ($$t-\tau$$). Depending on the values of the parameters, this equation displays a range of periodic and chaotic dynamics. Here, the rationale of this equation from a biological perspective, its mathematical properties, the equation's development, and some open problems are explored.

## Feedback systems

In feedback systems, the value of a control variable is sensed and appropriate changes are made in the rates of production and/or decay, with a typical goal of achieving a stable unchanging output. Physiological control systems represent prototypical feedback systems. Relevant variables that are subject to feedback control include the concentrations of electrolytes, oxygen, glucose, and blood cells in the blood, blood pressure to the brain and various organs.

Perhaps the simplest differential equation representing a feedback system is

$\tag{2} \frac{dx}{dt}= \lambda - \gamma x,$

where $$\lambda$$ and $$\gamma$$ are positive constants and $$x$$ is a variable of interest. Equation (2) represents a system in which $$x$$ is produced at a constant rate $$\lambda$$ and destroyed at a rate $$\gamma x\ .$$ For this equation, once an initial condition is specified at some time (typically $$t = 0$$), the solution is completely determined and in the limit $$t \to \infty\ ,$$ $$x \to \lambda/\gamma$$ with an exponential approach to the asymptotic value. Thus, there is a stable fixed point in this equation which is attained from any initial condition.

In real feedback systems, there is typically a time lag between the sensing of the value of a variable under control, and the mounting of an appropriate response. In physiological systems, the time lag can be substantial. For example, following a loss of blood cells, it can take many days before new blood cells can be produced through the activation, differentiation, and proliferation of the appropriate blood stem cells.

Thus, in more realistic models of physiological feedback systems, the rates $$\lambda$$ and $$\gamma$$ at time $$t\ ,$$ are not necessarily constant, but may depend on the value of $$x$$ at time $$t - \tau\ ,$$ designated $$x_{\tau}\ ,$$ where $$\tau$$ represents the time lag inherent in the control system. One of these equations, Equation (1), which displays a sequence of period-doubling bifurcations and chaotic dynamics, has become known as the Mackey-Glass equation (1977). In this equation, $$x$$ is a concentration (assumed non-negative for all times) of circulating blood cells, and $$\beta_0$$ and $$n$$ are constants involved in the description of the dependence of the production of these cells as a function of $$x_{\tau}\ .$$ The physiological rationale for the nonlinear production control function is as follows. If $$0<x_{\tau} \ll 1$$ the individual would be very sick and unable to generate many blood cells. If $$1 \ll x_{\tau}\ ,$$ the person would have too many blood cells so that the production rate would once again be low. For some intermediate values of $$x_{\tau}$$ the production rate would be maximal. Figure 1 shows the dynamics of Equation (1) with $$\gamma=0.1\ ,$$ $$\beta=0.2\ ,$$ $$n=9.65\ ,$$ $$\tau=2$$ starting from an initial condition of $$x=0.5$$ for $$t<0\ .$$ The dynamics are chaotic. Figure 2 shows a time-delay embedding of the dynamics. Figure 2: Time delay embedding for the Mackey-Glass equation, (1), for $$\gamma=1\ ,$$ $$\beta=2\ ,$$ and $$\tau=2\ ,$$ $$n=9.65\ .$$

## Relevance to nonlinear dynamics

Equation (1) is a simple equation that can generate complex dynamics including chaos. However, in contrast to low dimensional dynamical systems, such as the Lorenz equation (1963) and the Rössler equation (1979), time- delay differential equations such as Equation (1) are infinite dimensional systems. This is because it is necessary to specify an initial function over the time interval $$[-\tau,0]$$ in order for the solution to be well defined and to be able to integrate the equation. Doyne Farmer, who worked on the characterization of chaotic dynamics and was the first to call Equation (1) the Mackey-Glass equation, recognized that increasing the value of $$\tau$$ increases the dimension of the attractor in chaotic systems (Farmer, 1982). This observation, and the simplicity of the equation, has led to the evolution of this equation into one of the standard models used to test algorithms for the quantitative characterization of chaotic dynamics. Such analyses focus, for example, on the determination of the dimension of the attractor, the Kolmogorov entropy, the Lyapunov exponent, and predictability (Farmer, Ott, and Yorke, 1983; Farmer and Sidorowich, 1987; Grassberger and Procaccia, 1983).

The Mackey-Glass equation has also had an impact on more rigorous mathematical studies of delay-differential equations. In 1977, at the time of the original publication, methods for analysis of some of the properties of delay differential equations, such as the existence of solutions and stability of equilibria and periodic solutions had already been developed (Hale, 1977). However, the existence of chaotic dynamics in delay-differential equations was unknown. Subsequent studies of delay differential equations with monotonic feedback have provided significant insight into the conditions needed for oscillation and properties of oscillations (Mallet-Paret and Nussbaum, 1989; Walther, 1995; Mallet-Paret and Sell, 1996). For delay differential equations with non-monotonic feedback, mathematical analysis has proven much more difficult. However, rigorous proofs for chaotic dynamics have been obtained for the differential delay equation $$\frac{dx}{dt}=g(x(t-1))$$ for special classes of the feedback function $$g$$ (Lani-Wayder and Walther, 1996). Further, although a proof of chaos in the Mackey-Glass equation has still not been found, there continue to be advances in understanding the properties of delay differential equations, such as Equation (1), that contain both exponential decay and nonmonotonic delayed feedback (Röst and Wu, 2007). The study of this equation remains a topic of vigorous research.

## Relevance to physiology

Originally, Mackey and Glass presented equations of the form of Equation (1) to illustrate the appearance of complex dynamics in physiological control systems by way of bifurcations in the dynamics. They suggested that many physiological disorders, called dynamical diseases, were characterized by changes in qualitative features of dynamics (Mackey and Glass 1977; Glass and Mackey 1979). The qualitative changes of physiological dynamics corresponded mathematically to bifurcations in the dynamics of the system. The bifurcations in the equation dynamics could be induced by changes in the parameters of the system — as might arise from disease or environmental factors, such as drugs — or changes in the structure of the system. Subsequently, a large number of researchers have analyzed bifurcations in mathematical models of physiological systems and associated these with the abnormal (pathological) dynamics of disease. Interesting examples arise in hematology, cardiology, neurology, and psychiatry (Bélair, Glass, an der Heiden, and Milton, 1995; Milton and Jung, 2003).

## Searching for chaos

An interesting sidelight is the way Mackey and Glass searched for the chaotic dynamics. They found that by using a time delay embedding of the equations, in which the value of $$x$$ was plotted as a function of $$x_{\tau}\ ,$$ the plot would develop kinks prior to the initiation of chaos. Figure 3-12 shows a sequence of plots of these time-delay embeddings for different values of $$n\ .$$ Although there was not sufficient space in the original Science article to show these plots, they were shown by Glass at the New York Academy of Sciences meeting conference on Chaos in November 1977 (Mackey and Glass, 1979). At this meeting David Ruelle suggested to Glass that it would be interesting to look at the return map on a section through the attractor as represented in the time-delay embedding to see if a map similar to the quadratic map would be found. The resulting plots, which have not been published previously, do show interesting structure, similar to that anticipated by Ruelle (Figures 13 and 14). To the best of our knowledge, there have been no systematic studies of return maps computed in this way and a careful numerical and mathematical analysis of the properties of these return maps still seems warranted.

About the same time, Doyne Farmer and collaborators published their studies in which derivative embeddings and time-delay embeddings were used to help characterize the geometry of chaotic dynamics (Farmer, 1982; Farmer, Ott, and Yorke, 1983; Packard, Crutchfield, Farmer, and Shaw, 1980).

## History and open problems

The development of this equation reflects a convergence of research interests and collaborations between Mackey and Glass. Mackey went to the McGill Department of Physiology in 1971 from the NIH, and at that time was working on problems in membrane ion transport. His interests quickly shifted, however, to an examination of the dynamics of cell differentiation and replication as a result of his work with John Combs. Specifically, he became interested in the dynamics of hematological disease in which the pathology was associated with pronounced oscillations in blood cell numbers.

In 1975, Glass also took a position in the Department of Physiology at McGill. Prior to coming to McGill, Glass had been working on nonlinear dynamical models of biochemical and gene control systems, first at the University of Chicago and then at the University of Rochester. Theoretical studies of biochemical feedback systems with spatially localized catalysts and coupling by diffusion displayed oscillations even in cases where the spatially homogeneous system did not show oscillations (Glass and Kauffman, 1972; Glass and Pérez, 1974). The diffusion-induced time delays led to the stabilization of oscillatory dynamics. This suggested the possibility of inserting time delays in systems of differential equations, and Glass carried out some unpublished simulations in which time delays were inserted into ordinary differential equation models of negative feedback systems.

During the summer of 1975, just around the time he was moving to Montreal, Glass attended a summer workshop at the Aspen Institute for Theoretical Physics. As recounted by Gleick (1987), participants at the workshop included Steven Smale and Mitchell Feigenbaum and there was a good deal of discussion about chaotic dynamics in the quadratic map. Figure 14: The return map computing the minimum value of $$x(t)$$ for $$x(t)<0.8$$ as a function of the preceding minimum value of $$x(t)$$ for $$x(t)<0.8$$ using the data in Figure 2.

After Glass moved to McGill, Mackey and Glass initiated discussions and collaborations that were focussed on determining whether chaotic dynamics could exist in physiological systems, and if they might be playing a role in some of the complex rhythms observed in physiology. Mackey had already started modeling control of blood cell production, and Glass initiated studies of respiratory control. In both circumstances time delays arise naturally due to physiological mechanisms, so the use of time delay differential equations was a natural way to proceed. Further, a recently published paper by May and Oster (1976) developed the notion, based on the observations in Li and Yorke (1975), that non-monotonic feedback in difference equations (such as the quadratic map) could lead to chaos in a model of an ecological system. It seemed possible that non-monotonic feedback in delay differential equations could also lead to chaos. This behavior was sought, and found, in the models for hematopoietic control. Further, just as in the quadratic map, Mackey and Glass found a sequence of period doubling bifurcations leading to chaos. In the sequence shown in Figures 3-12, Figures 4 and 5 correspond to period 2, Figure 6 corresponds to period 4, Figure 9 corresponds to period 3, and Figure 10 corresponds to period 6. An open problem is whether there is a complete correspondence between the well known sequences of bifurcations in the quadratic map and in Equation (1).

Science abounds with examples of parallel discoveries of phenomena by isolated individuals or groups at virtually the same time. The Mackey-Glass equation provides another example. In the Fall of 1977 Mackey met Andrzej Lasota, a Polish mathematician, who had been working with Dr. Maria Ważewska-Czyżewska (a hematologist) on problems involving blood cell pathologies. During that meeting Lasota told Mackey about a recent paper of his (Lasota, 1977) in which he had studied the solution behavior of the equation

$\tag{3} \frac{dx}{dt} = \beta x_\tau^n e^{-x_\tau}-\gamma x, \quad \gamma,\beta,n > 0,$

in which the nonlinearity had the same non-monotone character as in Equation (1). Lasota, too, had discovered similar behavior in his equation as Mackey and Glass had. Everyone was unaware of the work of the other. However, that 1977 meeting led to a long collaboration between Lasota and Mackey (see Mackey, 2007, for a recent recounting of this). But this was not the only coincidence. Perez, Malta, and Coutinho (1978) also studied differential delay equations but as a model for oscillations in laboratory blowfly populations. One of their models was

$\tag{4} \frac{dx}{dt} = (\beta_0 - \beta_1 x_\tau) x_\tau-\gamma x, \quad \gamma, \beta_0,\beta_1 > 0,$

and the solution behavior that they displayed in one figure clearly showed the numerical characteristics of chaotic dynamics. Thus, the bifurcations and chaotic dynamics that have been intensively studied and characterized in the Mackey-Glass equation, might have been pursued in the Lasota equation or the Perez-Malta-Coutinho equation but for some small chances of fate.

Acknowledgments: Thanks to Raluca Apostu and Thomas Quail who provided technical assistance with the figures and text. Jianhong Wu and Roger Nussbaum provided useful comments on an earlier version.