The study of a dynamical system and its bifurcations can be used to investigate neuronal behaviour. The pupil light reflex has an important neuronal component and it can be observed non-invasively in well documented and well controlled conditions. There is a vast amount of published work on the pupil light reflex, its anatomy and importance in clinical neuroscience, and on the mathematical tools, like Delay Differential Equations, used to model it. The exact characteristics of the neuronal component of the pupil light reflex are still not fully known and mathematical modelling of the nonlinearities involved can bring a better understanding of the data obtained from measurements of the pupil response to light.

A biologically plausible mathematical model of the pupil light reflex with Delay Differential Equations, extending a model by A. Longtin, John Milton et al. from McGill University, Canada, is introduced and analysed in "Nonlinear shunting model of the pupil light reflex", by Paul Bressloff, Cristina Wood and Peter Howarth, Proc. Roy. Soc. B vol. 263, 1996 and Spontaneous oscillations in a nonlinear delayed-feedback shunting model of the pupil light reflex (in PDF format), by P. C. Bressloff and C.V.Wood, Phys. Rev. E, vol. 58, 1998. A number of neuronal units are modelled by delayed differential equations as leaky integrators with shunting, meaning that the rate of change in the membrane potential depends on the difference between the current membrane potential and the equilibrium value. This way the mathematical model implements the logarithmic-like saturation of a neuron's response to very large input. The time delay in the equations is the delay between the moment when the light stimulus reaches the retina and the neuron starts firing in response to this stimulus. The scope of the model is to obtain maximum information about the neuronal population from measurements of the pupil size in response to light stimuli.

There is a subtle relation between the pupil area that is determined the iris muscle activity, and the retinal area responding to light, that determines the light flux. In normal conditions, the light stimulus has a spatial field larger than the eye and the responding area of the retina is determined by the pupil size (closed loop). There are also Maxwellian view experimental conditions, in which the light stimulus is a spot narrower than the smallest size that the pupil can reach, and the responding retinal area does not depend on the magnitude of the pupil area (open loop).

The parameters form a multi-dimensional continuous space.
As they vary, the stability of solutions can change, a process studied by bifurcation theory for dynamical systems.
An aspect of interest is the onset of oscillations,
both slowly oscillating solutions (passing through the equilibrium point at time intervals larger than the maximum time delay),
and rapidly oscillating solutions.
Slowly oscillating solutions can be used to describe the pupil cycling at pupil-edge light stimulation
and rapidly oscillating solutions can be used to investigate the pupillary hippus (fluctuations in pupil size at steady illumination).

The behaviour of the model is compared to data from measurements of the pupil light response
for a better understanding of the parameters and nonlinearities involved and their importance in medical diagnosis.
There is a large amount of noise and variation in the experimental data and numerical simulations are not sufficient
in understanding the qualitative behaviour of the model in a multi-dimensional parameter space.

The model considers at this stage only the parasympathetic branch of the pupil light reflex, concentrating on the neural component. The equations do not take into account interactions between neurons, that are not very well known. The model is like a jigsaw-puzzle: as more biologically plausible features are added (nonlinearities, more neuronal units, interactions between neurons, interactions between the paths corresponding to each eye, the inhibitory activity related to the sympathetic branch, a more detailed modelling of the retina and of the iris muscle, using Bayes probability theory to estimate parameters), the model fits better the experimental data and its parameters and nonlinearities can be better mapped to the real components of the reflex loop, widely studied in medicine and ergonomics.

## Description of the pupil response to light

Diagram of the human eye with a schematic enlargement of the retina, from webvision.med.utah.edu

The light enters the eye through the pupil and it is converted to nervous activity by the retina, a
complex photosensitive structure at the back of the eye. The light signal travels via the optic nerve from the eye
to the visual cortex. Some of the optic fibers go to a relatively small neuronal cluster for each eye,
called the Edinger-Westphal nucleus,
that sends signals to the pupil iris muscle, thus controlling the size of the pupil.
The pupil contracts when the retinal illumination increases and dilates when the retinal
illumination decreases or when the eye adapts to moderate ambient light.
This mechanism protects the retinal photoreceptors from bright lights and keeps on average a constant
light flux F. This light flux is given by all F_{i} received by the retinal ganglions, with
F_{i}=I_{i}A_{i}
(I_{i} is the light intensity received by,
and A_{i} the receptive area of retinal ganglion i). The pathway described above is related mainly to the parasympathetic component of the reflex.
There is also a very important sympathetic component, related to response to emotions, ambient temperature, pain, etc.
The parasympathetic and sympathetic activities define the autonomic nervous system, and this makes the study of the pupil
response to light very relevant to clinical neuroscience.
It is important, for medical research, to be able to isolate the behaviour related to the parasympathetic
and to the sympathetic system in the pupil light reflex. The noise and variation in the experimental data and
nonlinearities are making this a difficult task.

There are some important details concerning the neuronal component of the reflex, that are yet not fully understood, like the interactions between neurons, the contribution of the sympathetic branch, the mechanism of the noise called pupillary hippus (fluctuations in pupil size at steady illumination), adaptation, the response to dark pulses or to the switch between black and white patterns of equal overall illuminance (the pupil contracts), etc. The iris muscle and the neurons involved in the reflex have different response characteristics, and this is widely used in measurements. Neurons involved in the pupil light reflex have a firing rate that can reach 60 to 100 Hz, while the iris muscle cannot follow oscillations in the light signal of a frequency much higer than 4 to 6 Hz. This indicates that high frequency components of the pupil response to light might be related to the neuronal behaviour. There have been studies of the pupil response to sinusoidally modulated light, to beats, and of the pupillary noise. The neurons have also a different response time delay, typically 50 to 80 milliseconds, compared to a larger time delay associated to the iris muscle (in the region of 200 milliseconds). Measuring the pupil response time delay is thus used to better understand the neural component. One way of doing this is examining pupil oscillations created by pupil-edge stimulation: the light stimulus is a narrow spot of light directed to a point near the edge of the pupil. When the pupil contracts, it becomes too small and the light spot is covered, then the pupil dilates and gets to the size where the light spot can enter the eye, so the pupil will contract again, and so on. The delay in the pupil response to the light being on or off determines the period of these oscillations. This pupil cycling is important in medical diagnosis, since some nerve deffects influence the way the pupil oscillates, and the range of the pupil threshold area (determined by the position of the light spot) for which it oscillates.

## Oscillating solutions of Delay Differential Equations

The behaviour of Delay Differential Equations has been analysed in detail.
There are many relevant web links and
some classic reference books are
*Delay Equations. Functional, complex and nonlinear analysis*, by O. Diekmann,
S.A. van Gils, S. M. Verduyn-Lunel and H.-O. Walther, Springer-Verlag, 1995,
*Introduction to functional differential equations*, by J. K. Hale and
S. M. Verduyn-Lunel, Springer-Verlag, 1993,
*Differential-Difference Equations* by R. Bellman and K. L. Cooke,
Academic Press, London, 1963.

The system of n+1
nonlinear delay differential equations, that is model for the pupil reflex response to light
with n neuronal units modelled as leaky integrators with shunting,
represents a dynamical system.
A point in the phase space of this dynamical system is a function defined on the interval
[`t _{M} , 0`]
and with values in

`, where`

**R**^{n+1}`t`is the maximum of all time delay parameters. This phase space, with the supremum norm, is an infinite dimensional Banach space.

_{M}The existence and stability of periodic solutions for autonomic Delay Differential Equations can be
studied with methods like those used for Ordinary Differential Equations, modified to suit the structure of
the infinite dimensional phase space of the DDEs.
One method is to use a geometric approach similar to the Poincaré map for ODEs, and
define a map `A` given by the intersection of the solution to a convex closed subset `K`
of the phase space. A closed orbit representing a periodic solution will contain a fixed point of this map.
The problem is that in the infinite dimensional phase space, the closed convex subset `K` contains
the equilibrium point, that is a fixed point of the map. A nontrivial, nonconstant orbit
does not contain the equilibrium point. It is necessary to prove that the
origin is an ejective point and that the map has at least another fixed
point than the equilibrium point that is not ejective. The existence of this point implies the existence of
a nontrivial, nonconstant periodic solution. This solution is stable, as in the case of the
Poincaré map, if the map `A`
has an eigenvalue `v` with `v < 1`.

Another method, that analyses the behaviour of solutions in a neighbourhood of the equilibrium point under certain conditions of continuity and differentiability, is using centre manifold theory and normal form reduction. The characteristic equation of the linearized system has eigenvalues the roots of a transcendental polynomial. There is an infinite number of complex roots with negative real part. The number of eigenvalues with positive real part is finite or zero and changes only when parameters vary such that an eigenvalue crosses the imaginary axis. Thus, the system can move from a region in the parameter space with stable equilibrium point (where all eigenvalues have negative real part) to a region where the equilibrium point is unstable and a nonconstant periodic solution can appear (Hopf bifurcation). Centre manifold theory takes into account the nonlinearities involved and can determine if this periodic solution represents a stable or unstable limit cycle. The centre manifold contains all the bounded and periodic solutions in a neighbourhood of a non-hyperbolic equilibrium point (that has pure imaginary eigenvalues) and is invariant under the semiflow of the system of equations. The centre manifold is tangent at the equilibrium point to the centre subspace generated by the pure imaginary eigenvalues. The centre subspace contains all the periodic solutions of the linearized system of equations and the centre manifold represents the distortion caused by nonlinearities of this centre subspace.

It is quite difficult to study global behaviour of solutions, especially if the system has multiple equilibrium points. The system of equations considered is quite complicated, with many parameters, time delays and unknown exact expression for nonlinearities. Some of these nonlinearities can be non-monotonous and with hysterezis, making the global behaviour of solutions very difficult to analyze, like the real biological process that is being modelled.