Derivation of a Floquet Formalism within a Natural Framework

Many biological systems experience a periodic environment. Floquet theory is a mathematical tool to deal with such time periodic systems. It is not often applied in biology, because linkage between the mathematics and the biology is not available. To create this linkage, we derive the Floquet theory for natural systems. We construct a framework, where the rotation of the Earth is causing the periodicity. Within this framework the angular momentum operator is introduced to describe the Earth’s rotation. The Fourier operators and the Fourier states are defined to link the rotation to the biological system. Using these operators, the biological system can be transformed into a rotating frame in which the environment becomes static. In this rotating frame the Floquet solution can be derived. Two examples demonstrate how to apply this natural framework.


Introduction
Many biological systems are influenced by an environment, which is periodic in nature, e.g. due to the amount of light in the circadian cycle, or the seasonal weather conditions. To describe the evolution of such biological systems the periodic forcing of the environment should be taken into account. From a mathematical perspective, the Floquet formalism provides tools to deal with such recurring patterns (Farkas 1994).
The Floquet theory deals with systems of linear differential equations with periodic coefficients, and can be used to determine the stability of equilibria (Floquet 1883). Application to the field of epidemiology has been proposed previously (Heesterbeek and Roberts 1995a, b), but the use of the Floquet theory in ecology and epidemiology remains limited today. If applied, it makes use of numerical solutions of the linearized system of periodic ODE's (Klausmeier 2008). This is a straightforward and easily implemented method, but does not provide analytical solutions for the evolution of the biological system. For this reason the gap between the description of the evolution of the biological systems and the Floquet theory remains.
To bridge this gap, we wish to translate the Floquet formalism into a natural framework: the periodicity is dictated by the nature of the environment and is an integral part of the description of the biological system. Thus, we design a natural framework for the description of periodic causes. Furthermore, we derive the Floquet solutions of the linearized system of periodic describing the biological system within this framework and define the Floquet ratio R T with the same threshold property as the basic reproduction number R 0 . We provide a recipe and apply it to two examples, showing that this algorithm is relatively simple and can easily be used in calculations to determine stability of a specific system.

Conceptual Description
In this section we describe the concepts of the natural framework and the recipe for application. The full analytical description is presented in the Sect. 3, but for anyone who is only interested in practical application, it should be sufficient to read this section and continue reading in the Sect. 4.

Key Elements of the Natural Framework for the Effect of Periodic Fluctuations
We are interested in the stability of the null equilibrium of the population, which is stable if the long term growth rate of the population is negative. Due to the rotation of the Earth, the environment of the population is forced to be periodic in time and subsequently the population growth is time dependent. The time dependency of the change in the population growth is in which vðu 0 ; tÞ is the population state, K the periodic growth operator reflecting the influence of the environment on the population and t the time since start of invasion, u 0 the initial phase, which depends on the observation moment in the cycle. Because in the context of this paper the Earth rotation forces the environment of the biological to be periodic in time, u 0 should be interpreted as a parameter value for a variable u describing the rotation. In our description we undo this premature substitution and make the rotation explicit. Consequently, the environment is described by the growth operator Kðu; tÞ: The periodicity of this operator is forced by the rotation of the Earth with the frequency x r = 2p/T r , in which T r is the known length of the period. The different frequency components of the environment could be distinguished by means of the Fourier expansion, and therefore we can expand this operator in a Fourier series (Oppenheim et al. 1983) Kðu; tÞ ¼ in which K l are the Fourier components of K and the Fourier operators F l are defined by (Boender et al. 1998) These Fourier operators are used to describe the modulation of the environment by the Earth rotation. We describe the rotation of the Earth in physical terms by an angular momentum operator L z , which is an infinitesimal rotation in the plane perpendicular to the axis of rotation (z-axis) (Alonso and Finn 1968).
The periodicity of the growth operator can be removed by means of a transformation to the rotating frame, rotating with the frequency x r (Boender et al. 1998). The Floquet growth operator, i.e. the static growth operator in the rotating frame, is The Floquet growth operator is time independent, with a correction for the rotation -i x r L z . The analytical derivation of this equation will be performed in the Sect. 3. The dominant eigenvalue k max of K F determines the growth. If the real part r max of dominant eigenvalue k max is larger than 0, long term growth of the population will occur. When the real part r max is smaller than 0, the size of the population will fluctuate but still fades away to extinction. This threshold is equivalent to a description using the dominant Floquet multiplier exp r max ð Þ (Heesterbeek and Roberts 1995a, b). For exp r max ð Þ\1 the population will go extinct, but for exp r max ð Þ[ 1 the population will grow. If the interval between the generations (T c ) is exact and time independent, it is possible to construct a reproduction ratio R ¼ exp r max T c ð Þ (Wallinga and Lipsitch 2007). This well-known threshold quantity separates population growth (R [ 1) and decline (R \ 1). However, the generation interval T c is often periodic as well, due to for instance temperature dependent mortality rates. Therefore we propose here a general threshold quantity, which is straightforward to determine: the Floquet ratio which is the average number of offspring of a single individual after one period. This threshold quantity separates population growth (R T [ 1) and decline (R T \ 1) in the same way as R.

A Recipe for the Calculation of the Floquet Ratio
To calculate the eigenvalues of K F , we wish to construct a Floquet matrix of the Floquet growth operator. To this end, we have to formulate the matrices of L z and Derivation of a Floquet Formalism 305 F l . L z measures the Fourier number n of the function exp inu ð Þ and therefore the matrix can be formulated as F l increases the Fourier number n of the function exp inu ð Þ to n ? l. Therefore the matrix for F 1 can be formulated as and the matrix for F l can be found by multiplication of this matrix l times, because F l = (F 1 ) l . According to Eqs. 5, 7 and 8 the Floquet matrix can be written as in which the identity matrix I describes the rotational term i x r L z . In Box 1 a recipe is presented how to apply this description. In the Sect. 4 we will describe two examples.
In the Sect. 3 we will present a comprehensive derivation of the theory. Readers who wish to focus on the application can continue reading in Sect. 4.
Box 1 Recipe to determine the Floquet ratio (1) Establish a time dependent growth matrix Kðu 0 ; tÞ (see Eq. 1). This matrix describes the growth of the system (e.g. the increase of population size or the transmission and recovery causing a change in the number of infected hosts) (2) Determine the Fourier components K -n to K n of the growth matrix K (see Eq. 2) (3) Formulate the Floquet matrix (see Eq. 9) (4) Calculate the eigenvalues of the Floquet matrix, if possible analytically or otherwise using a numerical approximation (5) Compose the Floquet ratio R T from the real part of the largest eigenvalue 3 Theoretical Derivation

The Natural Framework
The variable u, as introduced in previous section, is essential for the constitution of the natural framework, because this is the variable for the description of the Earth's rotation. For an arbitrary function f ðuÞ holds f u ð Þ ¼ f u þ 2pk ð Þwith k an integer. This function could be expanded in terms of the Fourier states similar to a Fourier expansion (Oppenheim et al. 1983).
From Eq. 10 it follows that the Fourier states are normalized, because for the special case f u ð Þ ¼ U n u ð Þ, the coefficients f l equal the Delta function d nl , which is 1 for n = l and 0 otherwise. Therefore, U n u ð Þ constitute an orthonormal base for the description of functions and operators. The unit impulse or Dirac Delta is one the functions we applied further on (Oppenheim et al. 1983). This function is defined by its transformation property of an arbitrary function g Z p By means of a expansion, this Dirac Delta function can be expressed in terms of the Fourier states (Boender et al. 1998) according to Eq. 10. Furthermore, the Angular momentum operator L z and the Fourier operator F p can be formulated with respect to the Fourier states. The eigenvalue equation for the angular momentum operator becomes (Alonso and Finn 1968;Boender et al. 1998) Therefore, U l ðuÞ are the eigenstates of L z . When F p are applied to U l ðuÞ, this leads to Derivation of a Floquet Formalism 307 This operator causes an increase in the order of the function, i.e. l is transformed into l ? p. We will use this property in the derivations to calculate the multiplication of Fourier components.
To finish the picture, we derive the relation of L z and F p , using the expansion of the exponential function. This leads to the following relation for the Fourier operator in which n is an arbitrary state. In this derivation we use the expansion of an exponential function, the property exp ix r L z t ð Þexp Àix r L z t ð Þ¼1 and the following relation This natural framework facilitates the description of the population dynamics. The population could be expressed in terms of the Fourier states U l u ð Þ. This population is modulated by the environment, which is expressed in terms of the Fourier operators F p . The environment is forced to be periodic by the rotation of the Earth, which is described by the angular momentum operator L z .

Derivation of the Floquet Solution within the Natural Framework
We apply this natural framework to determine the growth of a population in a time periodic environment. The Floquet formalism unravels the solutions into two components: the long term behavior of the system and separated this from the short term (periodic) behavior of the system (Floquet 1883). The method is straight forward, in the sense that the long term behavior of a system is described directly by the Floquet exponents (Farkas 1994). If the largest of the Floquet exponents is larger than 0 the population will grow. According to Eq. 1, The biological system is described by a population state vðu 0 ; tÞ and the effect of the environment on this biological system by a growth operator Kðu 0 ; tÞ, in which the initial phase u 0 is a constant. We transform the nondifferentiable constant u 0 into a differentiable variable u to enable a description of the effect of the rotation on the biological system. This phase variable u can be transformed back to the initial phase u 0 , accordingly in which Eqs. 11 and 12 are applied. The mathematical solution of Eq. 1 is presented by Floquet (1883). This solution is applied in physics to solve the Schrödinger equation for a period Hamiltonian operator, which describes the energy of a system leading to the Floquet Theory (Boender et al. 1998;Shirley 1965). We will now apply this theory to the field of population biology. We present the algebraic derivation of disentanglement of periodic changes and the long term growth.
At the start of the invasion one typical individual is present. This initial state will be denoted by v 0 . An invasion can start at any time during the season. Therefore, this state v 0 is independent of t and u 0 . According to Eq. 17, Eq.
in which the evolution operator Uðu; tÞ describes the development of the biological system in time, and is defined by In the following we elaborate on the operator part of Eq. 18, resulting in disentanglement of the growth operator Kðu; 0Þ and the rotating period x r . Substitution of the outcome in Eq. 18 transforms the whole equation into a rotation frame, in which the growth operator becomes time independent. As a consequence this equation can be solved. After substitution of this solution in Eq. 17, the variable u is transformed back into a constant u 0 , leading to an explicite expression for the population state vðu 0 ; tÞ.
We use the Eqs. 15 and 2 to transform Eq. 18 into a constant growth operator in a rotating frame, Thus, the growth operator Kðu; 0Þ is disentangled from the rotation with x r around the axis of rotation (z-axis), which is described by the exponential operator exp ix r L z t ð Þ. To transform the total equation into the rotating frame we substitute the Floquet evolution operator U F ðu; tÞ ¼ exp Àix r L z t ð Þ Uðu; tÞ into Eq. 20, leading to Because these expressions equal zero according to Eq. 18, this leads to a directly solvable equation in which K F is the Floquet growth operator, according to Eq. 5. Equation 21 provides a solution for the Floquet evolution operator, U F ðu; tÞ ¼ exp K F ðuÞt ð Þ . Using the transformed Eq. 21, the solution of Eq. 18 is given by This is the operator presentation of the Floquet solution within the natural framework. Periodic changes and the long term growth are disentangled in this solution, in which ix r L z describes the periodic changes and K F ðuÞ the long term growth.

The Threshold Quantity for Long Term Growth
Because the Floquet growth operator, K F , contains only the Angular momentum operator L z and the Fourier operator F p (see Eq. 22), the eigenvalues k j and eigenstates w j ðuÞ of K F can be calculated. To do so, the initial state should be expressed in terms of these eigenstates: In Eq. 10 it is shown that the eigenstates w i ðuÞ can be expanded in terms of Fourier states, f il U l ðuÞ, where f jl is the lth Fourier coefficient of the Fourier expansion for the jth eigenstate. Substitution of Eq. 24 in Eq. 17 gives vðu 0 ; tÞ which is analogue to the Shirley formula (Shirley 1965). In the derivation the eigenvalue equation of the angular momentum operator in Eq. 13 and the orthonormality of the Fourier states according to Eq. 10 are applied. Equation 1 is now given in terms of Fourier coefficients, f jl , and Floquet exponents, k j . If the real part r max of dominant eigenvalue k max is larger than 0, long term growth of the population will occur. When the real part r max is smaller than 0, the size of the population will fluctuate but still fades away to extinction.

Approximation Methods for the Calculation of the Growth Rate
To calculate the eigenvalues of the Floquet operator, we applied approximation methods for the eigenvalues of this operator, which are are regularly applied in physics (Leskes et al. 2010). In this section we follow this approach and apply this to biological systems. To determine the eigenvalues, we have to find the diagonalized Floquet operator D F , which has the following property in which D F the diagonalization Floquet operator K F is the diagonalized Floquet growth operator, and K 0 it's zero Fourier component. The diagonalization Floquet operator can be defined by Derivation of a Floquet Formalism 311 This can be used to expand Eq. 26 according to in which the commutator of two operators A en B: [A, B] = AB -BA is applied. From Eqs. 3 and 16 the following commutation relations can be deduced: in which is used F m F n = F m?n . In the diagonalisation procedure, we expand the operator S F = S F (1) ? S F (2) … into several orders. To diagonalize the Floquet growth operator in first order, the (off-diagonal) terms of Eq. 5 should be removed, leading to This removes the off-diagonal blocks of K F as a result of the relation according to Eq. 29. Substitution of this in Eq. 28 gives The two latter terms are the higher order off-diagonal terms. If [K p ,K n ] = 0 for arbitrary p and n, these terms would disappear and K 0 ¼ K 0 . If the size of K n is much smaller than the rotation frequency, i.e. K n j j ( nx r , then the higher order terms in Eq. 32 are negligible. If this is not the case, the diagonalization procedure should be extended to the second order leading to in which the same procedure is applied as was done in the first step in Eq. 30. Substitution of Eqs. 30 and 33 into Eq. 28 would lead to a Floquet growth operator which is diagonal up to second order. This procedure should be repeated until desired convergence is obtained. In the next section we will provide two examples where we apply this formalism for the analytical determination and the numerical determination of R T .

Application
In this section we illustrate the recipe for the determination of the Floquet ratio (see Box 1) by means of two simple models. For a more complex model we refer to a study of Rift Valley fever (Fischer et al. in prep). In this subsection we describe an ecological example: the establishment of an exotic species. The establishment of an exotic species requires that the equilibrium without the new species is unstable. Seasonality will have an effect on a new species. Especially seasonal temperature changes strongly affects population growth of ectothermic animals such as arthropods, fish, amphibians and reptiles. To illustrate our Floquet method, we will use the example of a species with a temperature dependent intrinsic growth rate, gðu 0 ; tÞ (see e.g. Delatte et al. 2009). The intrinsic growth rate is the difference between the birth rate and the mortality rate. For a model of this biological system Eq. 1 has the form, The average daily temperature follows a sinus function with the periodicity x r ¼ 2p T r , in which T r = 1 year. For this simple example, we assume that the intrinsic rate of increase is linear with temperature. The intrinsic growth rate of increase has the form gðu 0 ; tÞ ¼ g mean þ g amp sinðx r tÞ, where g mean is the mean growth rate and g amp the amplitude of the intrinsic growth rate in time.
4.1.2 Determine the Fourier Components K -n to K n of the Growth Matrix Kðu 0 ; tÞ The expression gðu 0 ; tÞ could be expanded in Fourier series À 1 2i g amp expðÀðx r tÞÞ þ g mean þ 1 2i g amp expðx r tÞ. As a consequence, the Fourier coefficients are K À1 ¼ À 1 2i g amp ; K 0 ¼ g mean , and K 1 ¼ 1 2i g amp , according to Eq. 2.

Formulate the Floquet Matrix (see Eq. 9)
In this case the matrix K n equals the number K n and the identity matrix I equals 1.

Calculate the Eigenvalues of the Floquet Matrix
Because K -1 , K 0 , and K 1 are coefficients in this example, they all commute, i.e. [K p , K n ] = K p K n -K n K p = 0. According to Eq. 20, the eigenvalues of the diagonalized matrix is Derivation of a Floquet Formalism 313 4.1.5 Compose the Floquet Ratio R T from the Real Part of the Largest Eigenvalue r max = g mean is the real part of the dominant eigenvalue. Therefore the threshold quantity could be constructed, according to 6 and the Floquet ratio becomes If R T [ 1 i, the null equilibrium is unstable. Furthermore, the population cannot grow if R T \ 1. Obviously for this simple example, the Floquet ratio is larger than one if the growth rate is positive.
in which x and y are the number of infected vectors and hosts, and l x (t) is vector mortality and l y is host mortality. In this example we investigate a seasonal effect with the periodicity x r ¼ 2p T r , in which T r = 1 year. The transmission rate is a function of the seasons by average daily temperature (24 h) due to maturation of the eggs, the size of the population peaks twice during the period, and each day has a fluctuation in vector activity.
in which b mean is the mean transmission rate and the amplitudes of the transmission rate are due to the temperature b temp and due to the size of the population b pop , and the daily amplitude of the vector activity is b day . The mortality rate is only a simple function of temperature with a mean l mean and amplitude l temp : 4.2.2 Determine the Fourier Components K -n to K n of the Growth Matrix Kðu 0 ; tÞ The function of the mortality rate can easily be expanded in a finite Fourier series using the relation sin x r t ð Þ ¼ exp ix r t ð ÞÀexp Àix r t ð Þ 2i (see previous example). The Fourier components for the transmission rate are ±1, ±2 and ±365. We neglect the daily rhythm of transmission (±365), because this will most likely have a small amplitude and is effective on a different time scale. Using Eqs. 24, 37 and 39 the different Fourier components K l can be calculated: The Fourier coefficients are now represented as matrices, according to Eq. 40.

Calculate the Eigenvalues of the Floquet Matrix
Diagonalization of this Floquet matrix will provide the eigenvalues of this matrix. In this example, the calculation of the eigenvalues needs to be approximated by numerical techniques. In principle the Floquet matrix is infinite. For numerical calculation, the matrix must be truncated to determine the eigenvalues. This truncation will introduce numerical errors. According to Eq. 26, eigenvalues are not unique numbers, but should follow the pattern k 0i n x r (Farkas 1994;Shirley 1965). If the truncated Floquet matrix is extended, higher orders of n should appear. Eigenvalues which do not follow this pattern, are not the sought values, but numerical errors caused by the truncation. The eigenvalues of the Floquet matrix are in general complex, because in the evolution of interacting subpopulations, oscillations could occur. The imaginary part of k 0 will reflect these oscillations. However, the real part of the dominant eigenvalue determines whether growth occurs. For this reason we focus on the real part. The maximum value of the real parts of the eigenvalue r max determines the evolution of the epidemic.

Compose the Floquet Ratio R T from the Real Part of the Largest Eigenvalue
The Floquet ratio is given in Eq. 36. If R T [ 1 i, major outbreaks could occur and if R T \ 1 only minor outbreaks could occur.
Derivation of a Floquet Formalism 315

Discussion
We derived a methodology to determine the stability of the null-equilibrium for biological systems which are periodically forced by the environment. This formalism can be used for simple models in which one or a few parameters are forced to change periodically in time and for complex models including many forced time-periodic parameters. For simple models analytical results for the eigenvalues of the Floquet growth matrix K F can be achieved. The real part of the dominant eigenvalues r max of K F leads straightforward to the Floquet ratio R T , which has the well-known properties of a threshold quantity. For complex models a numerical calculation of r max is required, while all other steps in the determination of R T remain analytical. The formalism can be used for systems with a stable external forcing and a known periodicity. Examples are seasonality and the circadian rhythm which force the biological system by temperature and light conditions, and which have a regular period (i.e. of 365 days for seasons or 24 h for the circadian cycle).
Several methodologies were previously put forward to describe forced periodicities in biological systems. Analytical solutions for specific problems are proposed by Bacaer and Guernaoui (2006). However, in their approach for each new problem new solutions need to be constructed using the simplicity of the system. A general numerical approach using Floquet theory is put forward by Klausmeier (2008). This is not limited to specific problems. However, this technique does not provide an analytical framework for theoretical considerations. Several attempts for such analytical framework were discussed by Heesterbeek and Roberts (1995a, b). They investigated the suitability of the Floquet multiplier as a threshold quantity for epidemic growth. According to them, this quantity is probably the easiest threshold quantity to calculate numerically. They were not very satisfied with this quantity, because at that time an analytical Floquet framework was not available, in which analytical solutions and analytical approximations for this quantity could be derived.
We provided such Floquet framework in this paper. It is an analytical framework, including both analytical and numerical approximation methods. Within this framework the Floquet ratio is the natural threshold quantity. We offer a straightforward 5 step recipe, by which analytical solutions or analytical approximations of this threshold quantity could be derived. If numerical solutions are required, this is only needed in the fourth step for the calculation of the eigenvalue of the Floquet growth operator.
Within this natural framework we focused on the Earth's rotation as the single driving force of the periodicity. In the future, this can be extended in two ways. Firstly, a forced periodicity from another cause could also be described within this framework. The phase variable u would have a different origin and therefor a different meaning, but the mathematical formalism remains the same. Secondly, we present an example containing the seasonal rhythm and the circadian rhythm and we dealt with that by ignoring the relative fast circadian rhythm. A combination of several forced periodicities could be treated within a multi mode Floquet approach (Leskes et al. 2010). For each periodicity a phase variable should be defined where after the framework presented in this paper provide the tools to describe such a system. Although the procedure to formulate the multi mode Floquet matrix is straightforward, it is laborious to operate due to the large number of terms.