# Mathematical Phase Model of Neural Populations Interaction in Modulation of REM/NREM Sleep.

1 BackgroundOver the course of a sleep period, healthy humans experience a cyclical alternation of REM (rapid eye movement) and NREM (non-rapid eye movement) sleep. Each type of sleep has unique characteristics. NREM sleep is divided into four separate stages that follow progressively as sleep cycles advance, the last two being the stages of deep sleep.

REM sleep occurs in cycles of about 90-120 minutes across the night, and it accounts for up to 20-25% of total sleep time, although the proportion changes with age [18], and in pathological conditions such as Parkinson's disease [24]. This type of sleep dominates the latter half of the sleep cycle, especially the hours before waking, and the REM component of each sleep cycle typically increases as the night goes on. As the name suggests, it is associated with rapid, and apparently random, side-to-side movements of the closed eyes. REM has two components; the phasic component is a sympathetically driven state characterized by rapid eye movements, muscle twitches, and respiratory variability, while the tonic one is a parasympathetically driven state with no eye movements [11].

Although sleep is a restful time, it involves a complex activation of brain circuits as showed, since 1950s, through experiments using electroencephalography [19]. In late 1970's, Hobson and McCarley [13] explained REM/NREM alternation as a result of the antagonistic role played by two neuronal populations: FTG and LC neurons [17]. The locus coeruleus (LC) population, located in the posterior area of the rostral pons in the lateral floor of the fourth ventricle, is composed of mostly medium-size neurons, that produce norepinephrine. The activity of LC neurons is high during the NREM sleep and the wake. On the contrary, the gigantocellular tegmental field (FTG) population, a group of nerve cells in the pons, shows a concentration of discharges during the REM sleep. On the basis of these observations, McCarley and Hobson proposed a reciprocal interaction model for sleep cycle oscillation modulated by LC and FGT activity [13], [17]. The mathematical formulation of the model is based on a Lotka--Volterra system, in which the FGT cells as analogous to the prey population and the LC neurons represent the predator population. Recent advances in the neural anatomy and physiology involved in the regulation of the sleep have encouraged more detailed mathematical models.

Commonly, the neurons that stop, or significantly decrease, their firing rate during the REM sleep are called REM-off neurons, while neurons that significantly increase their activity during REM sleep are known as REM-on neurons [9]. In the first part of the present study, we have modeled the alternation of REM and NREM sleep starting from the classical Kuramoto model. The idea is to consider the two REM-on and REM-off neuronal populations, respectively promoting and inhibiting the REM sleep, as oscillators. Then, we have compared the synchronization of the classical Kuramoto system and the reaction--diffusion space time Landau--Ginzburg model. This comparison is the most innovative aspect of the paper.

2 Methods

We consider the interacting neurons as biological oscillators. Many types of physical, chemical and biological oscillators share an astonishing feature: they can be described by a single phase variable [theta] [2]. In the context of tonic spiking, the phase is usually taken to be the time since the last spike train, i.e. the time-series electrical signals recorded from individual neurons in the brain, essentially the action potentials generated by neurons [14].

2.1 Basic model

We base our analysis on the physical model of Kuramoto [1], that provides a system that can model synchronisation and desynchronisation in groups of coupled oscillators. The Kuramoto model considers a system of globally coupled oscillators modeled by the equation:

[mathematical expression not reproducible], (2.1)

where [[phi].sub.i] = [[phi].sub.i] (t) is the phase of oscillator i, [w.sub.i] is the natural frequency of oscillator i, N is the total number of oscillators in the system and K is a constant referred to as the coupling constant.

The derivation of the Kuramoto system in [15] is based on the complex Landau--Ginzburg equation (see equation (2.4.15) in [15])

[[partial derivative].sub.t][PSI] = A[PSI] + b[DELTA][PSI] - c[PSI] [[absolute value of [PSI]].sup.2], (2.2)

where A is a diagonal matrix with pure imaginary coefficients (1), b, c are complex coefficients while the term [PSI] [[absolute value of [PSI]].sup.2] is a self--interacting term. If the imaginary parts of b and c become very large, then the equation is very close to Schrodinger self- interacting system. As it was pointed out (p. 20, [15]) a chemical turbulence of a diffusion-induced type is possible only for regions intermediate between the two extreme cases, when Imb and Imc are very small or very large. If we have a system of N oscillators, then it is natural to assume [PSI] = ([[psi.sub.1],..., [[psi].sub.N]), such that are complex valued functions. If the oscillators have fixed amplitude then we can define a specific manifold, for example

[mathematical expression not reproducible],

such that the evolution flow associated to the problem (2.2) leaves M invariant. The additional restriction

[absolute value of [[psi].sub.j](t,x)] = 1, [for all]j = 1, ..., N, [for allt [member of] R, x [member of] [R.sup.n]

requires corresponding interpretation of the Laplace operator in the Landau --Ginzburg equation as Laplace--Beltrami operator generated by the fibre structure of the manifold M. Such type of models are studied intensively and a typical treatment can be found in [20]; so we have to replace (2) [DELTA][[psi].sub.j] by

[mathematical expression not reproducible]. (2.3)

Choosing the scale b, c ~ 0, Kb = [mu] with real [mu] we get the equation

[mathematical expression not reproducible]

and it is easy to see that the flow leaves M invariant. Moreover, the ansatz [[psi].sub.j] = r[e.sup.i[theta]j(t)], where r is constant leads to the classical Kuramoto system

[mathematical expression not reproducible].

The complex-valued sum of all phases

[mathematical expression not reproducible] (2.4)

describes the degree of synchronization in the network. The assumption r is a constant is not true, when [mu] is complex -valued. The in-phase synchronized state corresponds to r ~ 1, [PHI] being the population phase (the average phase). In contrast, the incoherent state, [[theta].sub.j] having different values randomly distributed on the unit circle, corresponds to r ~ 0, and we suppose the average phase to vanish.

Even in the real case, when [mu] [member of] R, the synchronization degree can be evaluated by the aid of the quantity

[mathematical expression not reproducible].

If all frequencies [w.sub.j] are equal, then the synchronization depends on the initial value D(0) and it is possible to check the exponential decay of the measure D(t) of the synchronization as t [right arrow] [infinity]. The Kuramoto model assumes the following:

* All oscillators in the system are globally coupled, i.e. each oscillator interacts with equal strength with all of other oscillators in the system.

* Individually, the oscillators are identical, except for possibly different natural frequencies [w.sub.i].

We start with the study of two coupled oscillators: we refer to two neurons, one REM-on with phase [theta], the other REM-off with phase [phi]. According to equation (2.1), if N = 2, their interaction can be described through the following phase model:

[mathematical expression not reproducible], (2.5)

where [w.sub.1], [w.sub.2] are the constants of frequency deviation from the free-running oscillation and A is linked to the constant K/2.

The highly nonsinusoidal nature of neuronal discharge activity suggests that a nonlinear interaction between REM-on and REM-off neurons is expected. So, we also study the general case

[mathematical expression not reproducible], (2.6)

where g is a nonlinear function such that there is A [member of] R for which g(u) ~ Au in a small neighborhood of u = 0, i.e. g is approximately linear in a small neighborhood of the origin.

As shown in Appendix 5.2, the classical theory of Lyapunov [12] is not sufficient to analyze the solutions of the systems (2.5) and (2.6).

Furthermore, since we look for the qualitative behavior of the solutions, by summing the equations of the two systems, we get

d([theta] + [phi])/dt = [w.sub.1] + [w.sub.2],

which implies the conservation law

[theta] + [w.sub.1] = ([w.sub.2] + [w.sub.1])t + k.

We consider the difference of the phases u = [theta] - [phi] and we obtain, from the system (2.5):

[??] = [w.sub.1] - [w.sub.2] + 2A sin u. (2.7)

The condition [w.sub.1] = [w.sub.2] is necessary, but not sufficient for the synchronization of the model. If we assume [w.sub.1] = [w.sub.2] and the initial conditions are given by [theta](0) = [[theta].sub.0] + [[epsilon].sub.1] and [phi](0) = [[phi].sub.0] + [[epsilon].sub.2] (where [[epsilon].sub.1] and [[epsilon].sub.2] are positive constants that represent small perturbations of the initial data), and [[theta].sub.0] ~ [[phi].sub.0], then u(0) ~ [[epsilon].sub.1] - [[epsilon].sub.2].

Similarly, by replacing u in (2.6), we obtain:

[??] = 2g(u) = 2Au + O([u.sup.2]) (2.8)

with initial conditions u(0) = [[epsilon].sub.1] - [[epsilon].sub.2]. Condition (2.8) just means that [??] - 2Au is approximately quadratic in a neighborhood of the origin.

The qualitative behavior of the synchronization measure depends on the sign of the constant A. In Figure 1 we present numerical results for the case A = -1, 0.1 = [w.sub.1] [not equal to] [w.sub.2] = 0. The case A = -1, 0.1 = [w.sub.1] [not equal to] [[omega].sub.2] = 0 is represented in Figure 2, where the synchronization effect is manifested on a different scale.

The effect of changing the sign of A is given in Figure 3.

Other cases for the choice of the initial data parameters and A are represented in Figures 4 and 5.

We consider [theta] and [phi] as the average phases of the two populations or clusters of neurons. More precisely, we can suppose that all REM-on and REM-off oscillators are globally coupled each other. So, according to Kuramoto (2.4), we can associate to the REM-on population the average phase [theta] and to the REM-off population the average phase [phi]. It is reasonable to consider the case of a large number N of neuron clusters, but the numerical simulations and the arguments based on rigorous a priori estimates show that the total qualitative picture for N = 2 and N [greater than or equal to] 3 is quite similar. However, our main goal shall be the comparison between the synchronization of the classical Kuramoto ODE system and the reaction--diffusion space time Landau--Ginzburg model that is determined in the next subsection.

2.2 Diffusive model

Turning back to Landau--Ginzburg equation (2.2), we use the idea suggested in [15] and choose an appropriate scale for models manifesting a chemical turbulence of a diffusion-induced type. We choose b = 1 for simplicity and assume the self--interacting effect are neglected; so with c ~ 0 we have the system

[mathematical expression not reproducible]. (2.9)

Using now the ansatz

[mathematical expression not reproducible],

we get a system of the form

[mathematical expression not reproducible].

Assuming [PSI] is in-phase synchronized state (3), we can simplify the above system to the following one

[mathematical expression not reproducible]. (2.10)

Our goal shall be the study of the synchronization measure

[mathematical expression not reproducible]

for the evolution flow associated to the system (2.10). This model seems to be a reasonable attempt to model the discharge activity of REM-on and REM-off neuron clusters localized in bounded domains.

Again we simplify the model (2.10) assuming N = 2; so we have only two oscillating phases [theta] = [theta] (t, x) and [phi] = [phi] (t, x), where the variable x localizes the position of the cluster spike, x is periodic and x [member of] (0, 2[pi]).

We use the principle of Kuramoto synchronization for the two clusters of REM-on and REM-off neurons and we take [theta] as the average phase of the cluster REM-on and [phi] as the average phase of the cluster REM-off. For simplicity we limit our study to the case of [theta] as the average phase of the whole REM-on population and [phi] as the average phase of the whole REM-off population.

By studying localized interactions in network models [25], more precisely networks of electronically coupled cells, we can arrive, recalling (2.15), to models involving diffusive effect, represented by a second order derivative. We deduce that the interaction between two different groups of neurons is described by the following system:

[mathematical expression not reproducible], (2.11)

where A is an appropriate constant. We generalize (2.11) to the system:

[mathematical expression not reproducible], (2.12)

where g is a particular nonlinear function and [theta] and [phi] the average phases of REM-on and REM-off populations, respectively.

In particular we focus attention on this general case and we assume that g(u) = Au + h(u), where A [member of] R and h is a small perturbation. More precisely:

[mathematical expression not reproducible]. (2.13)

So, we rewrite the system (2.12) as:

[mathematical expression not reproducible].

If we set [mathematical expression not reproducible], we obtain:

[mathematical expression not reproducible].

By summing the two equations of the system, we get

[mathematical expression not reproducible] (2,4)

and, by subtracting:

[mathematical expression not reproducible].

By rescalling [mathematical expression not reproducible] the equation becomes

[mathematical expression not reproducible]. (2.15)

We add the initial conditions:

[theta](0, x) = [[theta].sub.0](x), [phi](0,x) = [[phi].sub.0](x), x [member of] [1.sup.1]

In the next section we analyze the stability of the solution of the system (2.15) and of the equation

[mathematical expression not reproducible].

This simplified scalar equation manifests the main properties of the multi particle model for the case of small initial synchronization dispersion (4).

3 Results

3.1 Basic model

Let us rewrite the system (2.6)

[mathematical expression not reproducible] (3.1)

with initial conditions [theta](0) = [[theta].sub.0], [phi](0) = [[phi].sub.0], where g is a non linear, continuous function such that there exists a constant A < 0, satisfying g(u) ~ Au in a small neighborhood of the origin.

The solutions of this system of ordinary differential equations can be decomposed in three different parts: a constant vector, a linear in time mode and a dispersive part, i.e. modes decaying (exponentially) in time [12]. More precisely, the solutions of (3.1) are:

[mathematical expression not reproducible],

where v(t) decays exponentially, i.e. there is M [greater than or equal to] 0 such that [absolute value of v(t)] [less than or equal to] M[e.sup.At] for every positive t. We conclude that, since A < 0, the solution u(t) of (2.7) and (2.8) is asymptotically stable.

If, on the contrary, A > 0, even starting with small initial data, we obtain asymptotically a constant difference between the solutions [theta] and [phi] of (2.6). More precisely, we suppose that g is a C1 real function with B a constant such that

[mathematical expression not reproducible]. (3.2)

Then the solutions [theta] and [phi]] of (3.1) are global and satisfy

[mathematical expression not reproducible].

3.2 Diffusive model

The manifold

[mathematical expression not reproducible],

obviously is not invariant under the action of the Landau--Ginzburg evolution flow associated with (2.9). However we can find an analogue of invariant manifold for this infinite dimensional dynamical system. We shall suppose for simplicity that [DELTA] is the self--adjoint realization of the Laplace operator on [L.sup.2](0, 2[pi]) with periodicity condition

f(0) = f (2[pi]), f'(0) = f'(2[pi]).

Lemma 1. If the initial data

[[psi].sub.j](0, x) = [f.sub.j] (x) (3.3)

belong to the Sobolev space [H.sup.1](0, 2[pi]), j = 1, ..., N, then one can find a unique solution

[PSI](t, x) [member of] C((-[infinity], [infinity]); [H.sup.1]((0, 2[pi]); [R.sup.N]))

to (2.9) with initial data (3.3). If the initial data satisfy

[mathematical expression not reproducible],

then for any t [greater than or equal to] 0 we have

[mathematical expression not reproducible].

This result shows that the domain

[mathematical expression not reproducible].

is invariant under the action of the evolution flow (2.9). Therefore, we can introduce the synchronization measure

[mathematical expression not reproducible]

for the evolution flow associated to the system (2.10).

This observation suggests that the general Landau--Ginzburg diffusion model (2.9) manifests stronger synchronization compared with the classical Kuramoto model. This prediction is verified by corresponding numerical simulations for the approximated Kuramoto diffusion model (2.10).

Now we observe that the solution of the equation (2.14),

[mathematical expression not reproducible]

has bounded [L.sup.[infinity]] and [L.sup.2] norms. Indeed, the solution of the equation is given by

[mathematical expression not reproducible],

where [DELTA] = [[partial derivative].sup.2.sub.x] is the Laplace operator on (0, 2[pi]) with periodic boundary data, and [e.sup..[DELTA]t] is the Schrodinger semigroup generated by [DELTA]. The function [mathematical expression not reproducible] is bounded as a consequence of the maximum principle ([22] and [23]).

The existence of a global solution is a consequence of the fact that [e.sup.[DELTA]t] is a contraction in any Lebesgue space [L.sup.p](0, 2[pi]). In the case of A < 0, we obtain a priori bound for the system

[mathematical expression not reproducible]

under the assumption (2.13). We reduce the construction of the solution to the determination of a fixed point for the integral equation

[mathematical expression not reproducible], (3.4)

where G is the differential operator -[DELTA] - 2A that generates the semigroup {[e.sup.-Gt], t [greater than or equal to] 0}. The idea to prove the existence of a fixed point for (3.4) is based on the application of the exponential decay of the semigroup [e.sup.-Gt] in the case of A < 0 [16]. As a consequence the solution u(t,x) of (3.4) decays exponentially with respect to t and we can deduce the same result of the basic model.

In the case of A > 0, we show that, starting with small initial data, the dispersive mode is transformed into a form that tends exponentially to a fixed positive constant. More precisely, we consider g with the assumption (3.2); then the solution of the system

[mathematical expression not reproducible]

with 0 [less than or equal to] [u.sub.0](x) [less than or equal to] B for every x in [S.sup.1], satisfies

0 [less than or equal to] u(t,x) [less than or equal to] B, [for all](t,x) [member of] [0,T] x [S.sup.1].

The equation of the system above has the stationary solutions u = 0 and u = B and:

1. The solution u [equivalent to] 0 is unstable;

2. The solution u [equivalent to] B is asymptotically stable.

We have found again, for the case A > 0, the same result of the basic model.

The key element is the stronger dissipation effect generated by the dissipation term in (2.10), i.e. with the term [DELTA][[theta].sub.j]

[mathematical expression not reproducible]. (3.5)

Indeed, we choose N = 2 for simplicity in (3.5); so we have the system

[mathematical expression not reproducible], (3.6)

and then we take special Fourier initial data

[theta](0, x) = [(1 + k).sup.1/2] cos(kx), [phi](0, x) = 2[(1 + l).sup.1/2] cos(lx). (3.7)

If k = l = 0, then the solution to the corresponding Cauchy problem (3.6) (3.7) is a solution to the Kuramoto ODE system. If k, l [much greater than] 1, then

[mathematical expression not reproducible]

with

[mathematical expression not reproducible].

Then Figure 6 shows the 3D graphic of the density

D(t,x) = [[absolute value of [theta](t,x) - [phi](t, x)].sup.2]

of the synchronization dispersion

D(t) = [[integral].sup.2[pi].sub.0] D(t,x)dx.

The stronger decay of the solution is clearly manifested.

The synchronization dispersion D(t) is represented on Figure 7, where two cases are compared: the red line represents [10.sup.11] *[D.sub.k=0 l=0](t) (i.e. the case of solution to Kuramoto system), while the blue one represents the case [D.sub.k=5,l=6](t). It is clear that in the high frequency case the synchronization dispersion is much stronger than for the case of pure Kuramoto system.

In the discussion below, we give a biological interpretation of the solutions for the basic and the diffusive models.

4 Discussion

In the present study, we have described the interaction between two neuronal groups, the REM-on and REM-off populations, through a phase model, focusing the attention about their synchronization. We have assumed [theta] and [phi] as the average phases of the two populations, according to Kuramoto theory, and we have related the synchronization and the desynchronization (synchronization loss) to the alternation of REM and NREM sleep. Across a night sleep, REM and NREM episodes oscillate, and the activity of REM-on and REM-off populations is not synchronized. In other words, when the most of REM-on neurons are synchronized, the average phase of the population REM-on produces relevant effects and the average phase of REM-off population is minimal; in this case REM sleep dominates, according to the observation that REM-off neurons cease firing during REM episodes [21]. A symmetrical situation occurs in order to generate the NREM sleep.

Our numerical simulations for A = -1 and A = +1, represented in Figures 1 - 5, show that the qualitative behavior of the synchronization measure depends on the sign of the constant A. When A > 0, a desynchronized activity of the two populations has been observed, and we can consider B as the maximal phase difference between the two populations, reached when one group predominates on the other. We have a REM episode when we start from [u.sub.0] = [[theta].sub.0] - [[phi].sub.0] > 0 (i.e. REM-on group is more synchronized than REM-off), because [theta] - [phi] [right arrow] B and [theta] [much greater than] so we can assume asymptotically [theta] ~ B and [phi] ~ 0. On the contrary, if we start from [u.sub.0] = [[phi].sub.0] - [[theta].sub.0] > 0, a NREM episode occurs. A positive difference between [[theta].sub.0] and [[phi].sub.0], already in initial conditions (even if small), could suggest a discontinuity, so our model suffices to explain the cyclic oscillation of REM and NREM episodes during the night, but not the sleep onset. The explication of an initial distance between [[theta].sub.0] and [[phi].sub.0] probably involves some aspects of the human sleep, not considered in our model. For example, in the last three decades, the two-process model, proposed by Borbely and Daan [3], [6] posits that a homeostatic Process S interacts with a process controlled by the circadian pacemaker (Process C), in the regulation of the wakefulness/sleep cycle. Also a role of neuropeptides has been conjectured among the mechanisms that underlie the initiation, maintenance, and exit of sleep and wake states [8].

In the present study, we have supposed that REM-on and REM-off oscillators are globally coupled. According to Kuramoto theory, we have linked to the REM-on population the average phase [theta] and to the REM-off population the average phase [phi]. When we have considered a number of neuronal populations greater than two, unexpectedly numerical simulations showed that the total qualitative picture for N = 2 and N [greater than or equal to] 3 is quite similar. Nevertheless, the assumption that the two oscillators are globally coupled is biologically strong, if one considers the complexity of interaction between brain circuits and structures necessary to generate REM sleep. REM episodes result from complicate links of GABAergic, cholinergic, and aminergic neurons which control the activity of glutamatergic reticular formation neurons [4], [5], [7].

Our second model was a reasonable attempt to describe the discharge activity of REM-on and REM-off neuron populations localized in bounded domains. In this sense, this model is more realistic. Presumably, there are several different groups of neurons that are selectively "on" and "off" in REM sleep, so [theta] and [phi] have to be interpreted as the mean of two selected "on" and "off" clusters. In fact, our numerical simulations (Figure 6) showed a greater synchronization degree than the pure Kuramoto system, suggesting a more evident interaction if we consider only two homogeneous cluster of neurons. Experimental work is necessary to deepen the result that, in high frequency case, the synchronization loss is much stronger than for the classical Kuramoto model (Figure 7). One hypothesis could be that the second mathematical model better captures the presence of different frequencies in REM and NREM sleep, as recorded by electroencephalography [10].

5 Conclusion

The two mathematical models suffice to describe the alternation of REM and NREM sleep during the night. In particular, the second model is an innovative attempt to represent the synchronization and desynchronization, bounding the domain of two selected "REM-on" and "REM-off" neuronal groups. Considering that the reaction--diffusion space time Landau--Ginzburg model seems to be more precise in modelling the synchronization, it could be useful to analyze the hypersynchronization of neurons observed in epilepsy. Experimental data will be necessary to test the proposed models.

http://dx.doi.org/10.3846/13926292.2016.1247302

Acknowledgements

The third author is supported by University of Pisa, project no. PRA-2016-41 "Fenomeni singolari in problemi deterministici e stocastici ed applicazioni"; INDAM, GNAMPA--Gruppo Nazionale per l'Analisi Matematica, la Probabilita' e le loro Applicazioni and Institute of Mathematics and Informatics, Bulgarian Academy of Sciences and supported by Top Global University Project, Waseda University.

http://dx.doi.org/10.3846/13926292.2016.1247302

References

[1] J.A. Acebron, L.L. Bonilla, C.J.P. Vicente, F. Ritort and R. Spigler. The Kuramoto model: A simple paradigm for synchronization phenomena. Reviews of Modern Physics, 77(1):137-186, 2005. http://dx.doi.org/10.1103/RevModPhys.77.137.

[2] D. Bhowmik and M. Shanahan. How well do oscillator models capture the behaviour of biological neurons? Neural Networks (IJCNN), The 2012 International Joint Conference on, pp. 1-8, 2012. http://dx.doi.org/10.1109/ijcnn.2012.6252395.

[3] AA. Borbely, S. Daan, A. Wirz Justice and T. Deboer. The two-process model of sleep regulation: a reappraisal. Journal of Sleep Research, 2016. http://dx.doi.org/10.1111/jsr.12371.

[4] R.E. Brown, R. Basheer, J.T. McKenna, R.E. Strecker and R.W. McCarley. Control of sleep and wakefulness. The Journal of Physiology, 92(3):1087-1187, 2012. http://dx.doi.org/10.1152/physrev.00032.2011.

[5] R.E. Brown and J.T. McKenna. Turning a negative into a positive: Ascending GABAergic control of cortical activation and arousal. Frontiers in Neurology, 11(6), 2015. http://dx.doi.org/10.3389/fneur.2015.00135.

[6] S. Daan and C. Berde. Two coupled oscillators: simulations of the circadian pacemaker in mammalian activity rhythms. Journal of Theoretical Biology, 70(3):297-313, 1978. http://dx.doi.org/10.1016/0022-5193(78)90378-8.

[7] C.J. Van Dort, D.P. Zachs, J.D. Kenny, S. Zheng, R.R. Goldblum, N.A. Gelwan, D.M. Ramos, M.A. Nolan, K. Wang, F.J. Weng, Y. Lin, M.A. Wilson and E.N. Brown. Optogenetic activation of cholinergic neurons in the PPT or LDT induces REM sleep. Proceedings of the National Academy of Sciences of the United States of America, 112(2):584-589, 2015. http://dx.doi.org/10.1073/pnas.1423136112.

[8] J.J. Ftaigne, K.P. Grace, R.L. Horner and J. Peever. Mechanisms of REM sleep in health and disease. Current Opinion in Pulmonary Medicine, 20(6):527-532, 2014. http://dx.doi.org/10.1097/MCP.0000000000000103.

[9] P.M. Fuller, C.B. Saper and J. Lu. The pontine REM switch: past and present. The Journal of Physiology, 584:735-741, 2007. http://dx.doi.org/10.1113/jphysiol.2007.140160.

[10] D.W. Gross and J. Gotman. Correlation of high-frequency oscillations with the sleep-wake cycle and cognitive activity in humans. Neuroscience, 94(4):10051018, 1999. http://dx.doi.org/10.1016/S0306-4522(99)00343-7.

[11] C.D. Harris. Neurophysiology of sleep and wakefulness. Respir Care Clin N Am., 11(4):567-586, 2005.

[12] Ph. Hartman. Ordinary differential equations. Wiley, New York, 1964.

[13] J.A. Hobson, R.W. McCarley and P.W. Wyzinski. Sleep cycle oscillation: Reciprocal discharge by two brainstem neuronal groups. Science, 189(4196):55-58, 1975. http://dx.doi.org/10.1126/science.1094539.

[14] E.M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. The MIT Press, London, 2007.

[15] Y. Kuramoto. Chemical Oscillations, Waves and Turbolence. Springer-Verlag, New York, 1984. http://dx.doi.org/10.1007/978-3-642-69689-3.

[16] A. Lunardi. Analytic Semigroups and Optimal Regularity in Parabolic Problems. Birkhauser, 1995. http://dx.doi.org/10.1007/978-3-0348-0557-5.

[17] R.W. McCarley and J.A. Hobson. Neuronal excitability modulation over the sleep cycle: A structural and mathematical model. Science, 189(4196):58-60, 1975. http://dx.doi.org/10.1126/science.1135627.

[18] W. Moreas, R. Piovezan, D. Poyares, L.R. Bittencourt, R. Santos-Silva and S. Tufik. Effects of aging on sleep structure throughout adulthood: a population-based study. Sleep Medicine, 15(4):401-409, 2014. http://dx.doi.org/10.1016Zj.sleep.2013.11.791.

[19] G. Moruzzi. Active processes in the brain stem during sleep. Harvey Lect., 58:233-297, 1963.

[20] A. Nahmod, A. Stefanov and K. Uhlenbeck. On Schrodinger maps. Comm. Pure Appl. Mathematics, 56:0114-0151, 2003.

[21] D. Pal and B.N. Mallick. Neural mechanism of rapid eye movement sleep generation with reference to REM-OFF neurons in locus coeruleus. Indian J Med, Res, 125:721-739, 2007.

[22] A. Pazy. Semigroups of linear operators and applications to partial differential equations. Springer-Verlag, New-York, 1983.

[23] M.H. Protter and H.F. Weinberger. Maximum principles in differential equations. Springer, New York, 1984. http://dx.doi.org/10.1007/978-1-4612-5282-5.

[24] T.M. Rodrigues, A. Castro Caldas and J.J. Ferreira. Pharmacological interventions for daytime sleepiness and sleep disorders in Parkinson's disease: Systematic review and meta-analysis. Parkinsonism & Related Disorders, 2016. http://dx.doi.org/10.1016Zj.parkreldis.2016.03.002.

[25] I. Tsuda, H. Kujii, S. Tadokoro, T. Yasuoka and Y. Yamaguti. Chaotic itinerancy as a mechanism of irregular changes between synchronization and desynchronization in a neural network. Journal of Integrative Neuroscience, 3(2):159-182, 2004. http://dx.doi.org/10.1142/S021963520400049X.

Appendix

5.1 Lyapunov analysis

We consider again the system (2.6):

[mathematical expression not reproducible].

We decompose any solution [mathematical expression not reproducible], where ([[theta].sub.s],[[phi].sub.s]) solves the linear system

[mathematical expression not reproducible]

and ([[theta].sub.r], [[phi].sub.r]) solves the system

[mathematical expression not reproducible]. (5.1)

In order to get informations of the qualitative behavior of the solution of the system (5.1), we linearize the equation, obtaining

[mathematical expression not reproducible],

i.e.

[mathematical expression not reproducible].

The eigenvalues of this matrix are 0 and 2A, so we have directions of possible stability (instability) of the orbits, when A [??] 0. This analysis is insufficient because, when a Lyapunov exponent is zero, is not clear the solution behavior.

5.2 Differential equation and semigroup property

A semigroup is a family of linear operators [(T(t)).sub.t[greater than or equal to] 0] with the following properties:

T(0) = I, T(t)T(s) = T(t + s), t, s > 0. (5.2)

For example, the solution of the linear system

[mathematical expression not reproducible],

where A is a real n x n matrix, is the function u(t) = [e.sup.tA] [u.sub.0] with the matrix exponential [e.sup.tA] = [[summation].sup.[infinity].sub.n=0] [t.sup.n][A.sup.n]/n!.

More generally, if S : X [right arrow] X is a bounded linear operator on a Banach space X and u : [[t.sub.0],T] [right arrow] X, the equation

[mathematical expression not reproducible]

has the solution [mathematical expression not reproducible] (n times).

The partial differential equation

[mathematical expression not reproducible] (5.3)

can be also represented by

u(t,x) = ([e.sup.St]u)(x), t [member of] [[t.sub.0], T],

where [e.sup.St] is a suitably defined semigroup [16]. This is what we did in Section 4.2.

It is known that, if [laplace](X) is the set of the bounded operators defined on the Banach space X, then

[mathematical expression not reproducible],

where [omega] is the supremum of the real parts of the elements of the spectrum of S (e. g. eigenvalues).

The semigroup property (5.2) is the abstract representation of the uniqueness of the solution of Cauchy problems: namely if u(t,to,uo) is the solution of (5.3), then [22]

u(t,[t.sub.0],[u.sub.0]) = u(t, [tau], u([tau], [t.sub.0], [u.sub.0])), [for all][tau] [member of] [[t.sub.0], t].

Paolo Acquistapace (a), Anna P. Candeloro (a), Vladimir Georgiev (a,c) and Maria L. Manca (b)

(a) University of Pisa, Department of Mathematics Largo Bruno Pontecorvo 5, 56127, Pisa, Italy

(b) University of Pisa, Department of Clinical and Sperimental Medicine Via Roma 67, 56126 Pisa, Italy

(c) Faculty of Science and Engineering, Waseda University 3-4-1, Okubo, Shinjuku-ku, Tokyo 169-8555, Japan

E-mail(corresp.): candeloro@mail.dm.unipi.it

E-mail: acquistp@dm.unipi.it

E-mail: georgiev@dm.unipi.it

E-mail: l.manca@med.unipi.it

Received September 7, 2015; revised October 28, 2016; published online November 15, 2016

(1) we assume that [iw.sub.j] are the diagonal elements of this matrix with [w.sub.j] [member of] R

(2) to be more precise one can see Theorem 2.2 in the above cited book of Kuramoto, where modified Schrodinger map equation is defined by the aid of the operator in (2.3)

(3) recall that this means [r.sub.j] ~ 1

(4) see the precise statement in Lemma 1 below

Caption: Figure 1. Synchronization loss in the case A = -1, 0.1 [w.sub.1] [not equal to] [w.sub.2] = 0. Graphics of the measure D(t) of the synchronization. Initial data [theta](0) = 0.01, [phi](0) = 0.02.

Caption: Figure 2. Synchronization loss in the case A = -1, 0.1 = [w.sub.1] [not equal to] [w.sub.2] = 0. Graphics of the measure D(t) of the synchronization. Initial data [theta](0) = 0.03, [phi](0) = 0.02.

Caption: Figure 3. Similar synchronization loss in the case A = 1, 0.1 = [w.sub.1] [not equal to] [w.sub.2] = 0. Graphics of the measure D(t) of the synchronization. Initial data [theta](0) = 0.03, [phi](0) = 0.02.

Caption: Figure 4. Synchronization is sensitive with respect to the choice of the sign of A. Synchronization decays for A = -1, [w.sub.1] = [w.sub.2]. Graphics of the measure D(t) of the synchronization. Initial data [theta](0) = 0.03, [phi](0) = 0.02.

Caption: Figure 5. Synchronization is sensitive with respect to the choice of the sign of A. We have synchronization loss for A = 1, [w.sub.1] = [w.sub.2]. Graphics of the measure D(t) of the synchronization. Initial data [theta](0) = 0.03, [phi](0) = 0.02.

Caption: Figure 6. Graphics of the density D(t, x) of the synchronization dispersion D(t). We have synchronization combined with very strong dissipation for A = -1.

Caption: Figure 7. The synchronization dispersion D(t) for the diffusion--Kuramoto problem (3.6) with initial data (3.7). The red line represents D(t) multiplied by the factor [10.sup.11], when k = 0,1 = 0, while the blue one represents D(t) in the case k = 5,1 = 6.

Please Note: Illustration(s) are not available due to copyright restrictions.

Printer friendly Cite/link Email Feedback | |

Title Annotation: | rapid eye movement, non-rapid eye movement |
---|---|

Author: | Acquistapace, Paolo; Candeloro, Anna P.; Georgiev, Vladimir; Manca, Maria L. |

Publication: | Mathematical Modeling and Analysis |

Article Type: | Report |

Date: | Nov 1, 2016 |

Words: | 6010 |

Previous Article: | An Extension of the Product Integration Method to [L.sub.1] with Applications in Astrophysics. |

Next Article: | Fuzzy-Presic-Ciric Operators and Applications to Certain Nonlinear Differential Equations. |

Topics: |