Eusebius 1 Doedel

Nonlinear Numerics

The objectives and some basic methods of numerical bifurcation analysis are described.

Several computational examples are used to illustrate the power as well as the limitations of these techniques. Future directions of algorithmic and software development are also discussed.

I. Introduction

The principal objective of numerical bifurcation analysis, as opposed to numerical simulation, is to compute continuum families of solutions to well-defined operator equations. Such computational results give a deeper understanding of the solution behavior, stability, multiplicity, and bifurcations, and they often provide direct links to the underlying mathematical theories.

The basic numerical method is continuation, the language in which algorithms are typically expressed and analyzed is functional analysis, and the fundamental theoretical tool is often the Implicit Function Theorem (1FT). For ordinary differential equations (ODEs), the numerical algorithms have attained a high degree of reliability. For introductions see Rheinboldt [START_REF] Rheinboldt | Numerical Analysis of Parametrized Non/inear Equations[END_REF], Seydel [START_REF] Seydel | From Equilibrium to Chaos. Practical Bifurcation and Stability A nalysis[END_REF], and Doedel et al. [START_REF] Doedel | Numerical analysis and control of bifurcation problems: (I) Bifurcation in finite dimensions[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems: (11) Bifurcation in infinite dimensions[END_REF], and for a com- prehensive literature overview see Allgower and Georg [START_REF] Allgower | Numerical path following[END_REF].

Here we highlight some representative basic algorithms used in the numerical analysis ofnonlinear equations. The emphasis is on ODEs, for which there exist several software packages in which these techniques have been implemented. We describe computational results for four different problems, namely, a singular perturba tion problem, a problem with "bursting" phenomena (Plant's model), a problem with derivative discontinuities (Chua's circuit), and a problem with homoclinic bifurcations (coupled Josephsonjunctions). These examples are used to illustrate the power as well as the limitations of the numerical techniques.

We conclude by describing a selection of problems which we think merit fu ture work, since resulting algorithms and software would find many applications. Some current efforts in these directions are indicated.

Continuation

Regular solution points

Consider the equation G(x) = 0, where G:Rn 1 1 ---+ Rn. Generically, there will be curves (''branches'') of solutions. which can be effectively computed by pseudo-arclength continuation [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF]. Given a solution x 0 , say, and the direction of the solution branch, x 0 , a nearby solution x 1 is computed by solving (xl-xo)*x 0 = L1s, for x 1 • Typically Newton's method is used, and the next direction vector x 1 is then easily computed as a by-product. It is easy to see why and when this continuation technique works.

In a more general setting, suppose that G(x 0 ) = 0, where G:PA 1 ~PA 2 . Here 81 1 and 81 2 are appropriate Banach spaces. Suppose that x 0 is a regular solution point, i.e.

%(~) = Span{cp} , %(~*) = {0},
where G~* is the adjoint of G~. Consider again the extended system ~(x) _ { G(x), } (x -x 0 )*cp-s.

Theorem

Near a regular solution point x 0 , the equation ~(x) = 0 defines a unique solution curve x = x(s).

Proof This follows from the 1FT, since is nonsingular.

•

We note that, in some cases, numerical continuation computations can be rigorous; see, for example, Kearfott [START_REF] Kearfott | An interval step control for continuation methods[END_REF].

Simple singular points

Next suppose that G(x 0 ) = 0, where, again, G:PA 1 ~!!4 2 . Suppose that x 0 is a simple singular point, i.e. %(~) = Span{.pb.p 2 }, %(G~*) = Span{t/t}. Consider the equations .

-lG(x o+cl<fJt+*f2CfJ2+ y )+J-Lt/l = 0,] F(y,Jl,t.~>£2) = Y 41t -0, y*cp2 = 0.

We have F(O,O;O,O) = 0. Theorem Near (y,,u;£~>£ 2 ) = (0,0;0,0), the equation F(y,,u;£t.t. 2 ) = 0 defines a unique solution manifold with y(O,O) = 0, Jl(O,O) = 0.

Proof This also follows from the 1FT, since the Fn!chet derivative (Fy(y ,jl;t I ,f2)lF.u(y ,jl;t I ,£2) at (y,Jl;tb£ 2 ) = (0~0;0,0), given by is nonsingular.

•

Assuming F (i.e. G) to be sufficiently smooth, we can derive the Taylor expansion

1 ,f,*(/! m m .!,*(/! m m 1 ,/,*(/! m m ( ) _ _ o/ XX't' I '1' I 2 + o/ XX't' ) '1' 2 + _ o/ XX't' 2 '1' 2 2 + 0( 3) Jl f.t.f2 -2 1/J*l/1 fl l/l*t/1 flf2 2 1/l*t/1 f2 f '
and similarly for y(t: 1 ,£ 2 ). To lowest order approximation, the solutions of G(x) = 0 then correspond to the level curves

(1/J*G~x(/)t q;l)d + 2(1/J*G~xq>l (/)2)£1£2 + (1/J*G~x(/)2(/)2)£~ = 0.
This is the Algebraic Bifurcation Equation (ABE; see Keller ( 6)) from which, in the case of distinct real roots, two distinct directions x~xo+£1(/)1 +t:2q>2 of bifurcating branches can be computed.

The above procedure, basically the Lyapunov-Schmidt procedure [START_REF] Chow | Methods of Bifurcation Theory[END_REF], also applies to higher-order singular points, and is the first step in the reduction to normal form. Normal forms are important in understanding the types of solutions that can exist near a singular point. However, these solutions often exist only in a very small neighborhood; see, for example, Wittenberg and Holmes [START_REF] Wittenberg | [END_REF]. Normal forms are therefore of little use in determining the actual global behavior of a given system. For numerical purposes the importance of the above theorem lies in the ABE that it produces. Furthermore, the basic idea, namely that of "extended system", is at the heart of most techniques for regularizing singular computational problems. In fact, the pseudo-arclength continuation algorithm is the simplest instance of this idea.

Linear algebra

Extended systems for computing and continuing singular points often give matrices in Newton's method that have a special structure which can be exploited. A general presentation of the numerical linear algebra aspects of extended systems can be found in Govaerts [START_REF] Govaerts | Computation of singularities in large nonlinear systems[END_REF]. An example of a system with such a special structure is

(~ 0 A* 1/1*
where the n + 1 by n +m matrix A has rank n, and where A* is the transpose of A. This system arises in the continuation of folds (m = 1), simple singular points (m = 2), and more generally in the continuation of m-parameter optima of some objective functional on a solution manifold (3, 4).

We now illustrate a solution procedure for the above system that requires only LU- decomposition (11) of the matrix A. Pivoting should be used. We have

~¥(A)= Span{cp,-}7!," JV(A*) = Span{l/1} .
The first equation is Ax = f-zc. To be solvable, we must have

1/l*f '7 - - ~ -1/J*c.
The solution component x then has the form x = xP + I r= 1 a;q>;. Here xP, { q>;}r= 1 , and 1/1 can easily be computed once A has been LU-decomposed.

From the second equation it follows that q> t

C ( xP + LJ-1 ~<pj) = q>f g, for i = 1 , ... ,m, or m L (cptCcpj)a 1 = cptg-cptCxP, i= l, ... ,m, j = 1
which we solve for the { a;}7!, 1 • Then m x = xP+ I a;q>;. i= I Now find a particular solution yP to A*y = g-Cx, for which the LU-decomposition of A can be used. Finally we have y = yP + {31/J, where, by the third equation,

1/J*(yP+ /31/1) = h, from which u h -tjl*yp J.' - 1/1 * 1/1 .
This solution algorithm also establishes the nonsingularity of the full matrix, given that A has rank n, that cf/::!.l(A), and provided that the matrix Q, with entries Q;j cptCcp 1 , is nonsingular. These conditions are generically satisfied.

Ill. Numerical bifurcation analysis of ODEs

Here we describe some basic ODE continuation problems and we discuss a particularly suitable discretization technique, namely orthogonal collocation.

Boundary value problems

The ODE boundary value problem addressed in Doedel et al. [START_REF] Doedel | Numerical analysis and control of bifurcation problems: (I) Bifurcation in finite dimensions[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems: (11) Bifurcation in infinite dimensions[END_REF] is of the form u'(t) = f(u(t) ,A.), !E[O,l], u( •)/( •,• )ERn, A.ERn;., with boundary conditions b;(u(O),u(l ),A) = 0, i = 1 ,2, ... ,nh, and integral constraints I q,{u(t),A) dt = 0, i = I ,2, . . . ,n,, where we assume that n;. = nb + nq-n + 1, which generically leads to one-dimensional continua ("branches") of solutions. 3.2. Periodic solutions A practically important example of a system of the above form arises from the numerical continuation of periodic solutions. A numerical scheme for computing stable and unstable periodic solutions is given by the time-scaled ODE

u'(t) = Tf(u(t),A),
where T denotes the actual period, subject to the periodicity condition u(O) = u(l ), and a phase condition

J: u(t)*U'(t)dt = 0.
This integral constraint, where u is a nearby reference solution, is a necessary condition for D(a)=J~l l u(t+a)-u(t)ll 2 dt to be minimized over a. After appropriate discretization, solutions of the above equations can be continued using the pseudo-arclength continuation technique. Asymptotic stability and bifurcation information can be extracted from the linearized Newton systems.

Connecting orbits

Another ODE boundary value system arises from the numerical continuation of connecting orbits, i.e. orbits that connect fixed points of a vector field. The basic algorithm, which can be presented in various forms [START_REF] Lentini | Boundary value problems over semi-infinite intervals and their numerical solution[END_REF][START_REF] Beyn | The numerical computation of connecting orbits in dynamical systems[END_REF][START_REF] Friedman | Numerical computation and continuation of invariant manifolds connecting fixed points[END_REF], consists of continuation of solutions to the following set of equations: the time-scaled differential equations

u'(t) = Tf(u(t),A), u( • ), f( •,•)ERn, AERn",
the stationary point equations and the eigenvalue-eigenvecto r and normalization relations fu(Wt.A)Vt ; = J1 1 ;VJi, vt,.vli = 1, i = 1, ... ,n"

where, for simplicity, we assume that fu(w 0 ,tl) has n 0 distinct real positive eigenvalues J1o; with eigenvectors v 0 ;, and that fu(wb2) has n 1 distinct real negative eigenvalues J1 1 ; with eigenvectors vli. (In practice, it is better to work with an orthogonal basis of the corresponding invariant subspaces £ 0 and£" or their orthogonal complements.)

One can also add a phase condition L (f(u,}.) -f(U ,f..))*f.(U))f(U)) dt = 0. This integral constraint, with reference solution (u)), is a necessary condition for minimizing D(a)=fbiiu'(t+a)-u'( t)il 2 dt over er, also using the fact that u"(t) = fu(u,).)u'(t) = fu(u ,A.)f(u,A.). Finally, we require u(O) to lie in the tangent mani- fold £ 0 at "distance" fo from w 0 and, similarly, u(l) in £ 1 at distance £ 1 from w 1 :

no no u(O) = wo+Eo I Co;Vo;, I di = I,

i = I i = I n 1 n 1 u( 1) = w 1 + f 1 I c 1 ;V ]j, I cL = I. i =I i = I
These last equations can be used to eliminate w 0 and w 1 from the earlier equations. This leaves n coupled differential equations subject to ne= 2n + (n +I )(n 0 + n 1 ) + 3 con- straints. In addition to u(t)ER 11 we have scalar variables A.eRn'•, f 0 , f 1 ER, J1 0 ;,c 0 ,.eR, V 01 ERn, fori= I, ... ,n 0 , and J1. 1 ,-,c 11 ER, v 1 ,ERn, fori= l, ... ,n 1 • The total number of scalar variables is n 1 • = nJ. + (n + 2)(n 0 + n 1 ) + 2. Generically we need nl = ne-n. Thus the num- ber of parameters required for a single heteroclinic connection, i.e. the codimension of the orbit, is n;,-n-(n 0 + n 1 ) +I. A branch of orbits requires nA +I free parameters. The integration time T is "large" and fixed. It can be shown, under appropriate trans- versality conditions, that the above scheme is well-posed. A particular case is the heteroclinic connection of two saddle points in R 2 .

Here n = 2

and n 0 = n 1 = 1, so nA = I. A branch of such orbits requires two problem parameters. Another example is the connection of a saddle to a node in R 2 • Here n = 2, n 0 = I and n 1 = 2, so ni .. = 0. A branch of such connections requires one problem parameter. We also have the homoclinic orbit, for which w 0 = w 1 and n 0 +n 1 = n, so that n;. = 1. Often such orbits can also be computed as the limit of periodic orbits as the period T---+ oo .

Discretization

There are now several software packages for the numerical bifurcation analysis of ODEs (see Allgower and Georg (5) for a survey), some of which are more convenient to use than others. However, it is generally recognized that those based on global solution techniques (as opposed to "shooting") and orthogonal collocation with piecewise polynomial functions (as opposed to standard finite differences) perform best on difficult problems. Here "difficult" refers to boundary or interior layers, strong relaxation phenomena, sharp fronts and similar phenomena. A selection of such problems will be presented in Section 4.

Introduce a mesh {0 = t 0 < t 1 < ... < t N = 1 }, and let ~11 -tJ-tJ-1 for I < j < N. Define f!J>J: = {phe C[O,I] :p"ilr 1 .r 1 _dE&m}, where f!J>m is the space of polynomials of degree less than or equal to m. The orthogonal collocation method with piecewise polynomials (15) consists of finding phEfYJ'h' and A such that pl,(z 1 .J = f(ph(zJ.i),A), j = 1, ... ,N, i = 1, ... ,m, and such that Ph satisfies the boundary and integral conditions, including the pseudoarclength equation. In each interval [t 1 _ 1 ,t 1 ] the collocation points z 1 ,; are the transformed roots of the mth orthogonal polynomial ("Gauss points").

This method gives a globally accurate solution, which is especially accurate at the main mesh points t 1 and for the scalar variables A ("superconvergence"; see de Boor and Swartz ( 16)). Moreover, there is a known "optimal" mesh selection strategy that minimizes the number of mesh intervals [START_REF] Ascher | Numerical Solution of Boundar_v Value Problems for Ordinary Differential Equations[END_REF]. This safeguards, to some extent, against spurious solutions and, combined with continuation, makes it easier to compute difficult solutions. The collocation method also makes it easy to handle integral constraints.

Finally, the method leads to a sparse Jacobian matrix in Newton's method, which can be efficiently solved, allows the computation of Floquet multipliers at little extra expense [START_REF] Fairgrieve | Floquet multipliers[END_REF], and gives reasonable performance when implemented on a parallel computing system (18).

Collocation for PDEs

Piecewise polynomial collocation has been extended to certain classes of partial differential equations (PDEs), but its use there has not been as widespread as for ODEs.

Here we briefly describe a particular collocation method that may offer advantages in bifurcation software for systems of nonlinear elliptic PDEs. Consider

~u = f(u), xEO. c RN,
where ~ is the Laplace operator with appropriate boundary conditions and integral constraints. We suppress explicit indication of parameters. Suppose the domain is an N-cube, which has been recursively subdivided into 2M subrectangles ("finite elements"). To each finite element associate appropriate boundary matching points X;, fori= 1 , ... ,n, interior collocation points Z;, for i = I, ... , m, and a polynomial p(x) of "order" n +m. Require that, at the points X; on any common boundary, the values of neigboring polynomials and their normal derivatives match. Also require each polynomial to satisfy the collocation equations

Example

For the Laplace equation ~u

(x) = f(x) in R 2
, with a uniform mesh, one central matching point per element boundary, and one central collocation point per element, an equivalent "finite difference" approximation will have the form

4 vi = L a,.pj+ f3if(z), j =I
where u; and v; are the approximate solution and its normal derivative, respectively, at the matching point X;. As a basis, choose {I,x,y,x 2 ,y 2 }. Then, for each finite element, one obtains a discrete set of equations -------------------------------------------

~ 0 . 0 -0.5 -1.0 -1.5 -2.0 ,------------===~ -2.5•~--~----~--~----~--~----~--~---=~--~--~ 0.
v 1 = (2u 1 -u 2 -u 4 )/2h + hf(z)/2, v 2 = (2u2-Ut-u 3 )/2h + hf(z)/2, V 3 = (2u 3 -u 2 -u 4 )/2h+hf(z)j2, v 4 = (2u 4 -u,-u 3 )/2h + hf(z)/2.
This method is second-order accurate.

•

In its more general formulation, this PDE collocation method allows high-order accuracy and facilitates nested dissection, which reduces computational complexity [START_REF] Doedel | On the construction of discretizations of elliptic partial differential equations[END_REF].

IV. Computational examples

A singular perturbation problem

Consider the ODE boundary value problem (20)

t u"-uu'(u 2 -l) -u = 0,
with u(O) = 3/2, u(I) = y. Written as a system, this becomes u 2' = ).[ulu 2 (u~-I) + ul]/t , with u 1 (0) = 3/2, u 1 (1) = y. The solution behavior for varying y, which can be studied analytically, is illustrated in Fig. 1. There are six qualitatively different solutions for y in the range [-5/2,3/2]. The computations were done with E = 10-5 , using orthogonal collocation with up to 200 mesh intervals, 4 Gauss collocation points per mesh interval, and (obviously) adaptive mesh selection. The artificial parameter A has been introduced for homotopy purposes. We start with A = 0, for which there is an exact solution, and with a moderate value off, say, E = 1. The first numerical continuation takes A from 0 to 1. In the second continuation A= 1 remains fixed, while f decreases from 1 to I0-5 , say. The third continuation takes y from 3/2 to -5/2, with both E and A fixed. This final continuation is able to locate all different solution types that occur in the indicated range of y values. The calculations resolve the interior layers (in x), follow them as y changes, and determine they-ranges where the transitions occur. Although these calculations show the effectiveness of the continuation and discretization algorithms, it must be noted that the calculation with E = lo-s requires a very large amount of computer time. This is especially true in the range where an interior front moves with y, and where the pseudo-arclength step must be very small in order for the Newton iterations to converge. For larger£, say E = 10-3

, the calculations are relatively fast.

Plant's model of bursting nerve cells

A system is said to exhibit bursting when it changes back and forth between a quiescent state and a rapidly oscillating state. A five-variable ODE model for bursting in nerve cells developed by Plant [START_REF] Plant | Bifurcation and resonance in a model for bursting nerve cells[END_REF] is given by

V= [gisi(V)y 1 +gTxT][VI-V]+[gKxi+gpc(Kp+c) -'][VK-V]+gL[VL-V], XT= [sr(V)-xr]/(CrxT), XK = [sK(V)-xK]/(CtxK), Y1 = [zt(V) -YI]/((r_ l'f), C = p[KcxT(Vca-V)-c].
Detailed expressions for the various functions in this model and for parameter values can be found in Plant's paper. We let g 1 be the bifurcation parameter.

Figure 2 shows a portion of a branch of periodic solutions. The quantity Norm, which is used as a convenient solution measure, is defined as where T is the period of the oscillation. Over subintervals of the displayed range of 9~> the qualitative behavior of the periodic orbits is characterized by the number of spikes. Between these intervals are very narrow regions where one spike is added as g 1 increases. Figures 3 and4, with scaled-time variable t as horizontal axis and with V(t) as vertical axis, show representative bursting solutions carrying three and four spikes, respectively.

Figure 5 shows a blow-up of the bifurcation diagram around the second narrow region. To the left of this region, there is one spike; to the right, there are two spikes per period. Note that there is no branching process involved in the transition from one spike to two spikes. There are, however, cascades of period-doubling bifurcations. A

V 50.,.----------------------------------------------------~ 0.
-25.

-50. typical orbit on the first period-doubled branch is shown in Fig. 6. Note that the orbit alternates between one spike and two spikes.

The main branch terminates on the right in a homoclinic orbit. A very high period approximation to this homoclinic termination point is labeled 49 and shown in Fig. 7. It still carries four spikes, which are indistinguishable in the scaled-time plot due to the ----------------------------------------------------.

25.

0.

-2 5.

-50 . V 20 . .r-----------------------------------------------------.

0.

-10 .

-20.

-30.

-40 .

-so.

-60.

- On the other hand, these calculations demand much computation time and, in extreme cases, become impossible due to ill-conditioning of the boundary value formulation. (Initial value calculations are entirely useless in resolving the transition regions and the homoclinic behavior.) Although the use of more numerically stable linear solvers may somewhat extend the computation domain, it is likely that the calculations described here approach the computability limits. Two time-scale geometric singular perturbation analysis is likely to be the remaining (analytical) tool. In fact, such analyses can be found, for example, in Rinzel [START_REF] Rinzel | Bursting oscillations in an excitable membrane model[END_REF] and Terman [START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF]. See also Decroly and Goldbeter (24).

70.~----r---~----~----~----r---~-----r----~----~--~ 0.
V 50•~-----------------------------------------------------, 25

Chua 's circuit

Chua's circuit is one of the simplest electronic devices to exhibit complex behavior. A reference related to the discussion here is Khibnik et al. [START_REF] Khibnik | Continuation techniques and interactive software for bifurcation analysis of ODE's and iterated maps[END_REF]. The equations modeling the circuit are where

x' = a[y-h(x)], y' = x-y+z, z' = -11 v fJ. ' 1 h(x) = m 1 x+ 2(m 0 -mi){Ix+ 11-lx-11},
and where we take f3 = 14.3, m 0 = -1/7 and m 1 = 2/7. Note that h(x) is not a smooth function, and hence the solution to the equations may have non-smooth derivatives. However, for the orthogonal collocation method to attain its optimal accuracy, it is necessary that the solution be sufficiently smooth. Moreover, the typical adaptive mesh selection strategies will fail if the solution, or one of its lower-order derivatives, has discontinuities. Although one can be very precise about the necessary degree of continuity, it is in practice reasonable to require that solutions be c oc•. On the other hand, the adaptive mesh selection enables the resolution of "near-discontinuities", i.e. narrow regions in which the solution or some derivatives rapidly change. For these reasons, it is desirable that discontinuities be CXJ approximated. Thus we use the smooth approximation 2x lxl ~arc tan( Kx ), n which get better as K-HX). In the numerical calculations, we have used a very large value of K, namely K = 10 6 • Figure 8 shows a period-doubling sequence and Fig. 9 displays the last computed orbit on the second period-doubling branch. All three branches terminate in Shilnikov homoclinic orbits and, due to symmetry, these orbits occur in pairs. The displayed orbit, which has period T = 2000, is an approximation to one such homoclinic orbit.

It is represented in Fig. 9 using x and x ' as coordinates, which makes the rapid but smooth change of x' near x = -1 clearly visible. The accuracy and adaptive mesh selection of the orthogonal collocation method enables the resolution of this nearsingularity, which occurs six times along the displayed orbit. Moreover, the perioddoubling bifurcations remain detectable along the computed solution branches. However, the continuation process and the Floquet multiplier computation become less robust as K increases. ---------------------------~~------------------- ---~----

I ' -_,_ --" " I lll! --- ' \ ----- --a -------- ----------~llll ----------- ' / / ' ' -~- ---- 0.01~--------~------~--------r---===-~--------~------~ 6.

7.

8 .

9.

10.

11.

12.

ALPHA Fig. 8. A portion of the bifurcation diagram for Chua's circuit.

X-PRIME

3.~----------------------------------------------------~ -4•~--------~------~--------~--===-~--------~------~ -2.5 -2.0 -1.5 -1.0 -0 .5 0 . 0 0. 5 X PERIOD 2000
Fig. 9. An approximate homoclinic orbit (period 2000) in Chua's circuit.

Coupled Josephson junctions

Linear arrays of current biased Josephson junctions which are coupled through a shared load are used in various physical devices [START_REF] Hadley | Phase locking of Josephson-junction series arrays[END_REF]. Results on the general case of n junctions can be found in Aronson et al. (27). Here we consider the case of two junctions treated in Aronson et al. [START_REF] Aronson | A codimension-two point associated with coupled Josephson junctions[END_REF]. Let <pj be the phase difference across the gap in the jth junction, for j = I ,2. If the load is a capacitor then, after some algebra and appropriate normalization, the evolution of the <pj is governed by the equations where {3 is the capacitance of the individual junctions and I is a bias current. Assume that IE(O, I), and let a= arcsin I and p = n-arcsin I. Then there are four rest points for (<r>~><f>~,<p2,<p;), namely,

3 -(a,O,p,O), 2= (p,O,p,O), 1 -(a,O,a,O), 4 = (p,O,a,O).
Here the rest points have been visually ordered as they appear in the projected (cp"cp 2 ) phase plane; see Figs 10111213for examples. Rest point 1, which is a sink, is not shown in Figs I 0--13 as it is of no interest here. Rest point 2 is a hyperbolic point with 2dimensional stable and unstable manifolds, and rest points 3 and 4 are hyperbolic points with 3-dimensional stable and !-dimensional unstable manifolds. We use a subscript + to denote translation by 2n in the first coordinate and a superscript + to denote translation by 2n in the third coordinate. Thus, for example, 3+ =(a + 2n,O,p,O), 3+ =(a,O,p+2n,O), and 3! =(a+2n,O,p+2n,O).

A solution { cp 1 (t),<plt)} is a rotation if there exists a minimal T> 0 such that <r>lt + T) = cp;{t) + 2n, for j = I ,2. An in-phase rotation is one for which <p 1 = cp 2 ; see numerical continuation studies, show that the general rotation solution structure is very rich and difficult to describe. In order to provide at least a skeleton of this structure, we have computed loci of rotations of infinite period, as they are generally in some way the limit of finite period rotations. In particular, all in-phase rotations originate from symmetric heteroclinic o. so',-----,----------------------------------.-----------.--. I connections between the rest points 2 and 2! . Since the stable and unstable manifolds of 2 and its translate 2! are 2-dimensional, this 2-+2! connection is codimension-1, and there is a curve I = /({3) in parameter space, namely the leftmost curve in Fig. 14, along which these heteroclinic connections exist. The existence of a codimension-2 bifurcation point, at solution 1 in Fig. 14, on the curve I= 1({3) is proved in Aronson et al. [START_REF] Aronson | Coupled arrays of Josephson junctions and bifurcations of maps with SN symmetry[END_REF]. This is a homoclinic twist bifurcation and there is a branch of asymmetric 2-+2! connections, called "bellows", that emanates from it. Figure 11 illustrates such a bellow. A different type of bellow, namely a 4-+4! connection, occurs along the curve which contains solution 3. Another solution type is the discrete rotating wave or "'POM", for "pony on a merrygo-round". We find two branches of infinite-period POMs. One contains 4-+3+ -+4! connections; see, for example, solution 4 in Fig. 12. The other consists of 2-+2+ -+2! connections; namely the branch carrying solution 5. There are also two branches of socalled infinite period semi-rotors. One consists of 4-+4+ connections (the branch carrying solution 6), while the other consists of 2-+2+ connections (the branch with solution 7).

All branches, except the in-phase branch, appear to emanate from a single point P, namely solution 8 in the enlarged diagram detail in Fig. 15. This is apparently a codimension-2 point at which there is a 3-+2+ , and by symmetry, a 4-+2+ connection; see Fig. 13. Generically, these connections are codimension-2, since the unstable manifolds of 3 and 4 are both one-dimensional, while the stable manifold of 2 is 2-dimensional.

Unfolding results in Aronson et al. (28) imply the existence of a countably infinite number of branches of 2-+2! connections near P that spiral into P, and a countably infinite number of branches of 4-+4! connections near P that terminate in P, but that do not spiral. The numerical location of such additional branches has been very difficult, especially the additional spiralling branches. In fact, we were at first unable to obtain numerical evidence of these theoretically predicted solutions. It was not until the search was extended to very large period that we were able to detect another such 4-+4! connection. (We mostly used the high-period •'periodic" approximation, rather than the asymptotic boundary conditions discussed earlier.)

V. Software development directions

There is a need for further refinement of existing continuation algorithms and software for bifurcation analysis, and there is a need for their extension to new classes of problem. Below, a brief overview is presented of some problems and applications that especially merit further attention. The main categories considered are algorithms and software for singularities and manifolds, for some special problems and for largescale problems.

Interactive software

There is much current interest in computer user-interfaces and communication tools. Although it is important that algorithms be made as accessible as possible, there is a danger of creating a virtual world of empty shells. The emphasis should therefore remain on the algorithmic development within simple structures that allow updating, extension, and comparison of algorithms. In dynamical systems, there has been some progress in this direction. One notable achievement is the DsTool software (29), which includes some continuation capabilities. Another ongoing project, CONTENT [START_REF] Kuznetsov | a multiplatform continuation environment[END_REF], based on the earlier LOCBIF (25, 31) and on AUTO (32), is primarily continuation oriented.

Singularities

5.2.1. Problems with symmetry. Multiple bifurcations do not occur in generic systems with one parameter. In fact, even simple bifurcation does not arise. Thus, one could concentrate on codimension-1 bifurcations, namely folds, Hopf bifurcations, perioddoubling and torus bifurcations, and homoclinic orbits, use extended systems to follow these in two-parameters, detect codimension-2 bifurcations, etc. Ironically, most systems that arise in applications are not generic. However, one can introduce suitable unfolding parameters, much like the variable J1 in Section 2.2, to regularize the computations. The same can be done to deal with Hamiltonian systems (for an example, see Doedel et al. (33)). Nevertheless, it is often useful to deal directly with non-generic systems, such as symmetric systems, as one may want to exploit symmetries to reduce the computational complexity. The software SYMCON [START_REF] Gatermann | Symbolic exploitation of symmetry in numerical pathfollowing[END_REF] is an important development in this direction. 5.2.2. Numerical linear algebra. In a typical numerical continuation calculation, most computer time is taken by the setting up and the numerical solution of the linear systems that arise in Newton's method. As mentioned above, if the problem has symmetries then the cost of the linear algebra can often be reduced. More generally, if the Jacobian matrices have a special structure then this can often be exploited in the linear algebra. An example that covers the two-parameter continuation of folds and higher-order singularities was given in Section 2.3. As important as efficiency is the numerical stability of the solution algorithms, given the nature of limiting problems encountered in continuation calculations. Further work on algorithms, and their inclusion in numerical bifurcation software, remains important. For recent work, see Govaerts (10).

5.2.3. Homoclinic bifurcations. In Section 3.3, we briefly described a scheme for computing branches of connecting orbits. The method applies to both homoclinic and heteroclinic orbits. Connecting orbits are important in the study of wave solutions to nonlinear diffusive systems. They also serve as "organizing curves" ("skeletons") from which other solution behaviors can be inferred, as in the Josephson junction example in Section 4.4.

Several unfolding results for singular homoclinic orbits have appeared in the literature; see Kuznetsov [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] for an introduction. Numerical determination and continuation of codimension-2 homoclinic bifurcations has rapidly advanced in recent years. See, for example, Champneys and Kuznetsov (36) and references therein. There is prototype software for the detection of 16 different codimension-2 bifurcations [START_REF] Champneys | HomCont: An AUT086 driver for homoclinic bifurcation analysis[END_REF] which has been incorporated in AUTO (19). 5.2.4. Optimization. There is much current interest in controlling chaos, as evidenced by the many articles on the subject. Continuation techniques have played little or no role in these investigations. However, continuation techniques can be very effective in the stabilization of a desired solution type. One approach is the successive continuation technique, a strategy for the location and continuation of quadratic singular points of successively increasing codimension [START_REF] Doedel | Numerical analysis and control of bifurcation problems: (I) Bifurcation in finite dimensions[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems: (11) Bifurcation in infinite dimensions[END_REF]. It applies in particular to the continuation of quadratic extrema of an objective functional on a solution manifold. Although, in practice, it is not possible to guarantee that a globally optimal solution will be reached, the technique does determine the fate of a given optimal solution as an additional parameter varies. This provides geometric insight. especially if there are few control parameters. Applications include, for example, the search for solutions with some desired feature, the stabilization of periodic orbits by controlling critical Floquet multipliers, and the optimal placement or annihilation of an instability source (e.g. a Hopf bifurcation). Although the procedure is more time consuming than standard optimization techniques, its Jacobian matrices have the special structure considered in Section 2.3. Efficient implementation in a continuation package remains an unaccomplished goal. 5.3. Manifolds 5.3.1. Implicitly defined 2D manifolds. The calculation of solution curves and the continuation of curves of singular points of successively increasing codimension give much insight into the solution structure of nonlinear equations. In some applications, however, it is useful to calculate two-dimensional, or even higher-dimensional solution manifolds or singularity manifolds. See Allgower and Georg [START_REF] Allgower | Numerical path following[END_REF] for a survey of some basic numerical techniques. A global scheme was developed by Henderson (38), and a current software project (Pisces) for computing manifolds and their singularities is in progress (39). 5.3.2. Invariant manifolds. An important manifold type is the invariant manifold in dynamical systems. Simple examples are the stable and unstable manifolds of hyperbolic fixed points. Software for these has been developed by, for example, Kevrekidis and Jolly (40) and Nusse and Yorke (41). Less work has been done on the computationally more difficult case of the stable and unstable invariant manifolds of a periodic orbit. This problem can be reduced, via the Poincare map, to computing the corresponding manifolds of a fixed point. This suggests shooting, but its poor performance for unstable orbits makes a direct approach preferable.

Particularly challenging is the continuation of an invariant torus, a generalization of the periodic orbit. Although a general persistence theory has been known for some time (42), feasible computational schemes have only recently emerged; see, for example, Lorenz (43) and references therein. The special case of an invariant "circle" has been treated more extensively; see, for example, Broer et al. ( 44) and references therein.

5.4. Special ODEs 5.4.1 . Delay differential equations. Relatively little work has been done on the development of numerical continuation and bifurcation software for delay differential equations and, more generally, for functional differential equations. There is a software package XPPAUT (45) that can do simulation calculations and some numerical continuation, but it does not include, for example, the continuation of periodic solutions. Indeed, the latter problem is inherently infinite-dimensional. Certain numerical schemes were developed and tested many years ago (46, 47), but no generally distributed continuation software exists for this purpose. 5.4.2. Differential algebraic equations. Differential algebraic equations (DAEs) have been of some interest in numerical analysis in recent years. For the case of importance here, namely ODE boundary value problems (Section 3), an orthogonal collocation code COLDAE has been developed by Ascher and Spiteri (50). However, a general continuation and bifurcation code for DAEs does not exist. 5.4.3 . Systems with discontinuities. As illustrated in Section 4.3, one can approximate discontinuous systems by smooth systems to which standard discretization and continuation applies. However, the "near-discontinuities'', which can be resolved by a good adaptive mesh, require a very accurate discretization, worsen the conditioning of the linear solvers, and reduce the continuation stepsize. In extreme cases, this will limit the effectivenes of a numerical analysis, and one must deal directly with the discontinuities. One approach is to treat the location of discontinuities as unknowns and, after appropriate transformations of the independent variable, obtain an enlarged system on the unit interval with scalar unknowns corresponding to the unknown interfaces. However, the number of interfaces may be unknown beforehand and change during continuation, causing implementation difficulties. There are also fundamental theoretical difficulties in the bifurcation study of such systems. 5.5. Large problems 5.5.1. ODEs. As already mentioned, there is a need for further refinement of the linear algebra algorithms in numerical continuation. This is especially true for large systems. Optimal mesh selection may also need re-evaluation. It is not clear, whether the current selection criteria are indeed optimal near, for example, homoclinic orbits. For very large ODE systems, there is a need for easy-to-install and easy-to-use parallel software versions. A prototype implementation is described in Wang and Doedel (18). 5.5.2. PDEs. The greatest challenges lie, not surprisingly, in the development of numerical continuation and bifurcation software for PDEs. There is such a package for scalar nonlinear elliptic PDEs on general domains in R 2 , namely PLTMG (51). It includes continuation and some branch switching capabilities, and has been very successful in various applications. The underlying solution technique is that of multigrid, which is very fast, but whose performance typically degrades when the problem becomes difficult. However, such iterative techniques are a necessity for large scale PDE problems in R 3 and higher dimensions. Good results have also been obtained with stabilized simple iteration schemes for computing stationary PDE solutions "with mostly stable modes"; see, for example, Schroff and Keller [START_REF] Shroff | Stabilization of unstable procedures: The recursive projection method[END_REF].

There remains a need for general bifurcation software for systems of elliptic PDE, subject to general boundary conditions and integral constraints. For the case of such systems on simple domains in R 2 , the generalization of the collocation method of Section 3.5 carries some promise. To become comparable to current ODE bifurcation software, it will be necessary to use adaptive meshes. In this case, the direct solution of the linear systems arising in Newton's method remains feasible, so that a high degree of robustness is possible. Interesting applications include, for example, 2D Navier Stokes problems. Extension to mixed PDE/ODE systems would be desirable.

VI. Discussion

The effective use of numerical continuation for the bifurcation analysis of nonlinear equations demands a certain level of knowledge of bifurcation theory, of numerics and, of course, of the intended application area. In the case of dynamical systems, this is in contrast with simple simulations, which will give results, whether correct or not, under almost any conditions. In a numerical analysis, one must understand and confront the fundamental reasons of a qualitative change of behavior, have some understanding of their possible impact on the numerics, and be able to interpret the results in the context of the application. Computational difficulties encountered, or even the results themselves, may not be understood or appreciated. without some knowledge of the fundamentals of bifurcation theory and numerical analysis. It seems desirable, there- fore, that applied dynamical systems in curricula have a firm theoretical basis, in terms of bifurcation theory, and a sound algorithmic basis, in terms of numerical analysis.

As the examples in the previous sections have illustrated, there are practical computability limits, in terms of computer time, necessary precision, and complexity of the solution structure. Often one is interested in a limiting behavior, for which continuation is a suitable tool, but which leads to ever increasing accuracy requirements and numerical ill-conditioning, hence requiring ever increasing machine precision and computer time. Thus, the significance of computed results cannot be fully understood without recourse to analytical methods, such as singularity theory, perturbation analysis, or global geometric methods.

However, within the limits of computability, the advantages of a numerical bifurcation analysis can be numerous. Transitions from stable to unstable behavior can often be pinpointed, while phenomena that underlie transitions can frequently be determined, whether physically observable (i.e. asymptotically stable) or not. The numerical results often provide a much simplified representation of the fundamental features of solution behavior, a "skeleton", from which other behaviors <.;an, to some extent, be deduced. Sometimes the numerical results may even suggest new methods or theories, as was the case, for example, for the coupled Josephson junctions model in Section 4.4. For dynamical systems, continuation greatly reduces the dimensionality of the computation space, as the results are not functions of initial conditions. Indeed, a carefully computed bifurcation diagram often carries more fundamental information than an entire picture book of simulation results can provide.

Ul 1 . 5 ,

 15 _==-

Fig. 2 .

 2 Fig. 2. Part of a branch of periodic solutions for Plant's model.

3 -Fig. 3 .

 33 Fig. 3. V(t) versus scaled time t for solution 3.

4 -Fig. 4 .Fig. 5 .

 445 Fig. 4. V(t) versus scaled time t for solution 4.

  value of the actual period. Nevertheless, a projected phase diagram (not shown) would clearly indicate four spikes. This orbit is another illustration of the accuracy and robustness of the computational methods in resolving near-singular problems.

Fig. 7 .

 7 Fig. 7. V(t) versus scaled time t for solution 15.

PERIOD 15 .

 15 0•~-

Fig. I 0 Fig. 11 .Fig. 12 .

 01112 Fig. I 0 for an example. Although the physical interest is mostly in in-phase rotations, it is also of interest to classify other solution types. Numerical simulations, as well as

Fig. 13 .

 13 Fig. 13. Solution 8.

Fig. 14 .

 14 Fig. 14. Loci of connecting orbits for the Josephson junctions model.

BETA 0 .IFig. 15 .

 015 Fig. 15. A blow-up of a neighborhood of the codimension-2 point.

  Promising recent work has been done by Cometto et al. (48) and Luzyanina et al. (49).


Silvia Iacobelli
email: silvia.iacobelli@chu-reunion.fr

Francesco Bonsante

Pierre-Yves Robillard

Néonatologie

Comparison of risk factors and perinatal outcomes in early onset and late onset preeclampsia: A cohort based study in Reunion Island

come

Comparison of risk factors and perinatal outcomes in early onset and late onset preeclampsia: A cohort based study in Reunion Island

Introduction

It is well recognized that preeclampsia is a heterogeneous syndrome and that the clinical characteristics of the early-onset preeclampsia (EO-PE, gestational age at delivery < 34 weeks) and lateonset disease (LO-PE, ≥ 34 weeks) are different [START_REF] Sibai | Pre-eclampsia[END_REF]. The pathophysiology of the two disorders is also supposed to be different: Vatten et al. have suggested that preterm delivered preeclampsia associated with low birthweight may be caused by underlying placental abnormality, whereas pre-eclampsia delivered at term may represent a mixture of conditions, ranging from mild pre-eclampsia with moderate placental affection to hypertensive conditions in pregnancies without placental dysfunction. [START_REF] Vatten | Is pre-eclampsia more than one disease?[END_REF]. Recently, Redman et al. have suggested that the LO-PE form also relies on the placenta, but that the placenta perfusion problems at term are due to a different placenta problem than the dysfunctional spiral artery remodelling pathway, which is more typical of the EO-PE pathophysiology [START_REF] Redman | IFPA senior award lecture: making sense of pre-eclampsiatwo placental causes of preeclampsia?[END_REF][START_REF] Redman | Preeclampsia, biomarkers, syncytiotrophoblast stress, and placental capacity[END_REF].

Several publications have described the associated maternal morbidities, birth outcomes, clinical and laboratory features of the two preeclampsia types [START_REF] Ny | Comparison of indications of pregnancy termination and prognosis of mothers and neonates in early-and late-onset preeclampsia[END_REF][START_REF] Kucukgoz Gulec | Comparison of clinical and laboratory findings in early-and lateonset preeclampsia[END_REF][START_REF] Lisonkova | Incidence of preeclampsia: risk factors and outcomes associated with early-versus late-onset disease[END_REF]. The mother's younger age at first birth [START_REF] Robillard | Etiology of preeclampsia: maternal vascular predisposition and couple disease?mutual exclusion or complementarity?[END_REF] and some maternal risk factors (such as metabolic syndrome and hypertension) have been addressed for having important roles in development of the EO-and the LO-form, respectively. Also, maternal, foetal and neonatal complications are different in EO-PE compared to LO-PE [START_REF] Aksornphusitaphong | Risk factors of early and late onset preeclampsia[END_REF][START_REF] Lisonkova | Maternal morbidity associated with early-onset and late-onset preeclampsia[END_REF].

In general, within a population, both forms are represented, with the vast majority of women developing LO-PE (and so, mild neonatal signs) and a minority presenting EO-PE (and so, prematurity and impaired perinatal outcome) [START_REF] Chaiworapongsa | Pre-eclampsia part 1: current understanding of its pathophysiology[END_REF]. However, we recently showed that predominance of early-or late-onset preeclampsia may have huge geographical differences as the rate of EO-PE (defined as preeclampsia delivered before 34 weeks of gestation) was strikingly high (31.2% of all preeclampsia) in a cohort from Reunion Island, compared to other countries (8-10%) [START_REF] Iacobelli | Pre-eclampsia and preterm birth in Reunion Island: a 13 years cohort-based study[END_REF].

Reunion Island is a French department in the Southern Indian Ocean. The peculiarity of this tropical region lays in the multiethnic origin of inhabitants [Africa and intermixed population (50%), Europe (27%), India (20%) and China (3%)]. Compared to Europe and mainland France, there is a younger reproductive population (the mean age at primiparity is approximately 23 years). Finally, accessibility to maternity services is easy, and high-quality care is provided free of charge by the French healthcare system.

The aim of the present study was to compare the clinical differences, the maternal risk factors and the pregnancy outcomes between EO-PE and LO-PE in this particular population of preeclamptic women with a very high prevalence of preterm delivery.

Materials and methods

Design and study population

"Southern Reunion" identifies a specific geographic and health area of our department, with 5400 deliveries on average per year, occurring in three facilities (2 primary care and 1 tertiary care center). In this territory, all preeclampsia are referred to the tertiary care maternity of Saint Pierre University hospital.

This was an observational cohort study of all singleton preeclamptic pregnancies who delivered in this third level center over 15 years (January 2001 to December 2015).

Clinical data and variables

The study sample was drawn from the hospital perinatal database, which prospectively records data of all mother-infant pairs since 2001. Information is collected at the time of delivery and at the time of infant hospital discharge and regularly audited by appropriately trained staff; for the purpose of this study, records have been validated and have been used anonymously. Maternal characteristics and clinical risk factors examined for potential association with preeclampsia included: maternal age, parity and gravidity, marital status, education, smoking or alcohol assumption during pregnancy, infertility treatment, prepregnancy body mass index (BMI), diabetes mellitus, chronic hypertension prior to pregnancy, coagulopathy, hypercholesterolemia, history of preeclampsia, history of renal or thyroid disease and congenital malformation. The following delivery and pregnancy outcomes were analyzed: medically induced delivery, caesarean section, abnormal foetal heart rate, postpartum haemorrhage, intrauterine foetal death [(IUFD) in utero or intrapartum death of a foetus delivered at 22 weeks of gestation or later], medical termination of pregnancy (MTP), early neonatal death (death between 1 and 7 days of life) and perinatal death (defined as IUFD + MTP + early neonatal death), HELLP syndrome, and eclampsia. The neonatal outcomes of interest were: gestational age at birth, birth weight, sex, Apgar@1 min ≤ 3, small for gestational age, (SGA, according to [START_REF] Brenner | A standard of fetal growth for the United States of America[END_REF] and need for intubation at birth.

Definition of exposure

Preeclampsia was defined according to the international recommendations which were used during most of the period covered in this study: systolic blood pressure ≥140 mmHg or ≥diastolic 90 mmHg at bed rest on at least two occasions 6 h apart, and proteinuria (excretion of ≥0.3 g/day, or at least 1+ on a urine reagent strip), measured after the 20th week of pregnancy (No authors listed, 2000).

Preeclampsia outcome was categorized in two groups: EO-PE and LO-PE defined as delivery before 34.0 and ≥ 34.0 weeks, respectively.

Statistical analysis

Categorical variables were presented as frequencies and continuous variables as means ± SDs.

Bivariate comparisons were performed using χ2 test, or Fisher exact test for qualitative variables and ANOVA or Kruskal-Wallis test when appropriate.

First, women with preeclampsia were compared to women without preeclampsia. Second, EO-PE were compared to LO-PE. In a second stage of analysis, we estimated the effect of dichotomizing gestational age at delivery before 37.0 and ≥ at 37 weeks on maternal characteristics and risk factors. Comparative analysis of delivery and neonatal outcomes was adjusted for gestational age at birth when appropriate.

Because at bivariate analysis preeclamptic women had significantly higher rates of pre-pregnancy BMI and BMI > 30 kg/m 2 , and following the previous observation of the impressive increase of obesity in pregnant women at Reunion Island [START_REF] Robillard | Rapport 2001-2014 Du Relevé épidémiologique périnatal Sud Réunion[END_REF], in a supplementary analysis, we plotted the trend of obesity rate with that of both EO-PE and LO-PE over the study period in the observed population. A p value < 0.05 was considered significant.

Ethics

This study was conducted in accordance with French legislation. As per to the French law a specific approval of an ethics committee is not required for this non-interventional study based on retrospective, anonymized data and written consent is not needed.

Results

Epidemiology of preeclampsia

Over the study period 62,330 deliveries occurred in the maternity department and were recorded in the database. The incidence of preeclampsia in this entire cohort was 2.37%, with rates of 2.29% and 6.85% in singleton and multiple deliveries respectively.

For the purpose of the study, only singleton deliveries (N = 61,062) were analyzed. The total number of preeclampsia over the study period was 1397 and the frequency of EO-PE and LO-PE was respectively 0.75 (464 of 61,062) and 1.5 (933 of 61,062) per 100 deliveries. Among all preeclampsia, 33.21% were EO and 66.79 were LO. Among alive births from preeclamptic delivery, 60% were preterm (born less than 37 weeks of gestation) and 40% were at term (≥37 weeks of gestation). Among all infants 413 (31%) were < 34, 386 (29%) were 34-36 and 534 (40%) were ≥ 37 weeks of gestation.

Sociodemographic characteristic of the study population and comparison between EO-PE and LO-PE

Table 1 shows sociodemographic characteristics of the study population. Preeclamptic women were in general older and more often primiparae than no preeclamptic ones. The only difference when comparing EO-PE and LO-PE was the older age of primiparae with EO-PE.

When dichotomizing preeclampsia into preterm and term delivered groups (< 37 and ≥ 37 weeks of gestation, respectively) we found statistically significant differences between the groups for marital status and education (lower rates of "single" marital status and lower education level in EO-PE compared to LO-PE, p respectively = 0.03 and 0.04, data not shown). Sociodemographic characteristics of preeclamptic women were not different in patients who delivered a SGA infant (data not shown).

Maternal and clinical risk factors in the study population and comparison between EO-PE and LO-PE

All the considered gestational morbidities, as also obesity and higher pre-pregnancy BMI were significantly more prevalent in preeclamptic woman, compared to no preeclamptic ones, with no difference between EO-PE and LO-PE excepted for a rate of chronic hypertension significantly higher in EO-PE (Table 2).

When dichotomizing preeclampsia into preterm and term delivered groups (< 37 and ≥ 37 weeks of gestation, respectively) all the trends for women characteristics and risk factors remained similar and we found the following statistically significant differences between the two groups: higher rates of preexisting diabetes and lower rates of obesity in EO-PE (p respectively 0.04 and 0.02, data not shown).

Delivery outcomes in the study population and comparison between EO-PE and LO-PE

Delivery outcomes were different in preeclamptic women compared to no preeclamptic ones. The risk of medically induced delivery and that of caesarian section were higher in EO-PE compared to LO-PE, even after adjustment for gestational age (Table 3).

Foetal and neonatal outcomes in the study population and comparison between EO-PE and LO-PE

Gestational age and birth weight were lower in infants from preeclamptic mothers compared to no preeclamptic ones and the former presented significantly higher rates of SGA. Pregnancies complicated by EO-PE or LO-PE were more often associated with unfavorable foetal and neonatal outcomes. After adjustment for gestational age, foetal and neonatal outcomes did not differ significantly between EO-and LO-PE. However, the rate of MTP was significantly higher in EO-PE. (Table 4).

Trends of EO-PE and LO-PE over the study period and comparison with the obesity trend in the study population

As shown in Fig. 1, the rates of EO-PE and LO-PE did not change in our population during the study period, despite a significant parallel increase of the obesity rate.

Discussion

This study aimed to compare the clinical differences, the maternal risk factors and the pregnancy outcomes between EO-PE and LO-PE in Southern Reunion Island. Despite some specificities of this population, in particular the very high prevalence of preterm birth from preeclamptic mothers already reported [START_REF] Iacobelli | Pre-eclampsia and preterm birth in Reunion Island: a 13 years cohort-based study[END_REF] and confirmed here, the present study failed to identify any particular difference when comparing the two types of preeclampsia. The most impressive difference between the two groups refereed to the obvious more unfavorable neonatal outcome in infants born from EO-PE mothers, but this was merely dependent from the impact of a more preterm delivery.

Among socio-demographic factors, none conferred a relatively higher risk for EO-PE (as opposed to LO-PE), with the only exception of older age in primiparae. These findings are not consistent and even antithetic with previous works which emphasized that the early-onset disease in primiparae is prominent in adolescent and younger unmarried women [START_REF] Robillard | Etiology of preeclampsia: maternal vascular predisposition and couple disease?mutual exclusion or complementarity?[END_REF]. Other investigations [START_REF] Stubert | Clinical differences between early-and late-onset severe preeclampsia and analysis of predictors for perinatal outcome[END_REF][START_REF] Madazli | Comparison of clinical and perinatal outcomes in early-and late-onset preeclampsia[END_REF], having a sample size smaller than our, also failed to identify differences in maternal age when comparing EO-PE and LO-PE.

Several maternal risk factors of preeclampsia in our population were pre-existing morbidities that also represent risk factors for cardiovascular disorders (hypertension, diabetes, obesity, renal disease and hypercholesterolemia). All these diseases have long been considered as classical markers for preeclampsia at antenatal booking (Milne et al., 2005; [START_REF] Steegers | Pre-eclampsia[END_REF] and one recent review clearly underlined that metabolic syndrome and pre-existing hypertension are more common risk factors in women with late onset disease compared to those with the early-one [START_REF] Chaiworapongsa | Pre-eclampsia part 1: current understanding of its pathophysiology[END_REF]. On the contrary, in our population women with EO-PE were not exempted from these chronic systemic conditions and compared to LO-PE they even exhibited higher rates of pre-pregnancy hypertension. This finding is not new, and it is consistent with that of a population-based study from Lisonkova et al., 2014, which analyzed 14,201 preeclamptic women among 456,668 singleton deliveries in Washington State, as also with the results of a case-control analysis in preeclamptic women in Thailand [START_REF] Aksornphusitaphong | Risk factors of early and late onset preeclampsia[END_REF]., We found a statistically significant difference for chronic hypertension when comparing EO versus LO PE in our cohort (p = 0.02). This result may simply reveal the severity of the underlying disorder in preeclamptic mothers that, as a consequence, are more likely to have a preterm delivery. Mode of delivery was more often medically-induced and by caesarean section in preeclamptic pregnancies compared to their counterpart, with significantly higher rates in EO-PE women, and this association remained valid once adjustment for gestational age was made. In one previous study [START_REF] Iacobelli | Pre-eclampsia and preterm birth in Reunion Island: a 13 years cohort-based study[END_REF] we showed that indication to delivery was compliant to guidelines in 92% of all preeclampsia in our maternity care, with 44% and 48% of, respectively, maternal and foetal indications. So, the very high rates of delivery induction due to medical decision in our cohort may be the reflex of the preeclampsia severity in Reunion island population. Moreover, the prevalence of eclampsia was similar within the two groups.

In our cohort, and as reported by others [START_REF] Madazli | Comparison of clinical and perinatal outcomes in early-and late-onset preeclampsia[END_REF][START_REF] Elvedi Gasparovic | What affects the outcome of severe preeclampsia?[END_REF][START_REF] Stubert | Clinical differences between early-and late-onset severe preeclampsia and analysis of predictors for perinatal outcome[END_REF], EO-PE showed substantially higher risk for adverse foetal and neonatal outcomes, compared to LO-PE and to no preeclamptic population. However, the newborn prognosis in terms of early neonatal death was not affected once adjusted for gestational age at birth. The difference in prognosis for perinatal death seems so probably conditional on higher rates of MTP in EO-PE.

It worthwhile to note that among all preeclampsia in our population, 33% were EO-PE, thus giving birth before 34 weeks of gestation and that among LO-PE, 29% of women delivered between 34-36 weeks of gestation. As a consequence, even if neonatal outcome do not seem to be affected by preeclampsia per se, these results deserve concern, as preeclampsia accounts for a consistent number of preterm births in our population, which confers a significant health care burden. Indeed, infants born less than < 34 weeks of gestation, like late preterm infants (34 through 36 weeks of gestation), are at higher risk of morbidity and mortality than term infants.

Another interesting finding of our work was that rates of preeclampsia were stable over the study period in our population, despite the alarming increase of obesity in pregnant women in Southern Reunion Island. Obesity is nowadays a recognized risk factor for preeclampsia [START_REF] Jeyabalan | Epidemiology of preeclampsia: impact of obesity[END_REF][START_REF] Poorolajal | The association between body mass index and preeclampsia: a meta-analysis[END_REF]. However, the mechanism underlying the link between preeclampsia and obesity is not completely understood, and the lack of parallelism between these two trends is another example of the heterogeneity of preeclampsia, especially in terms of risk factors.

For many years preeclampsia has been considered a complex, heterogeneous syndrome, described as a two-stage disease with different pathophysiology and clinical symptoms. With increasing knowledge, it appears that such a clear distinction in two stages is somewhat inadequate [START_REF] Redman | The six stages of pre-eclampsia. Pregnancy Hypertension[END_REF]. From a pathophysiological point of view, for instance, it has recently been addressed that late-onset disease also is associated with findings of syncytiotrophoblast stress, which have been traditionally reported in the early-onset form [START_REF] Redman | Preeclampsia, biomarkers, syncytiotrophoblast stress, and placental capacity[END_REF]. From a clinical perspective, our study shows that several symptoms more often associated to the late-onset clinical syndrome (insulin resistance, perturbed coagulation function, hypercholestoremia) are also found in EO-PE.

As our study is a retrospective analysis, it has some notable limitations. Some interesting data are lacking in the perinatal database, such as the notion of primipaternity, length of sexual relationship, recurrent pre-eclampsia, maternal inflammatory state, thrombophilia and we acknowledge that this information could add an useful contribution to the exploration of risk factors associated to EO or LO subtype.

The strength of our study includes a large cohort sample from a specific geographic area which is very representative of the regional population.

Conclusion

In conclusion, even if compared to LO-PE EO-PE is a pathophysiologically distinct and more severe clinical entity with earlier gestational age onset and delivery, this study failed to identify a specific phenotype of preeclampsia in terms of predisposition or pre-existing risk factors for one form or another. Of course the complexity of the disorder and the wide range of interactions among risk factors associated to preeclampsia warrant further investigations. Few additional significant differences in terms of associated maternal morbidities appeared between the two groups when the analysis dichotomized for preterm versus term delivery, but these were not striking. Future research should be targeted on the identification of a more appropriate gestational age cut off for defining EO and LO-PE, with a special focus on factors influencing neonatal outcome in infants born from preeclamptic mothers.

Table 1

 1 Sociodemographic characteristics of the study population and comparison between EO-PE and LO-PE.

	Variables	No PE N =59,665 N (%) LO-PE N =933 N (%) EO-PE N =464 N (%) LO-PE vs No PE p	EO-PE vs No PE p value EO-PE vs LO-PE p
					value		value
	Age, years *	27.6 ± 6.6	28.5 ± 7.0	29.2 ± 6.8	< 0.0001	< 0.0001	ns
	Age, years * (primipara only) 23.7 ± 5,5	24.8 ± 6.0	26.1 ± 6.6	< 0.0001	< 0.0001	0.02
	Adolescent (12-17 years)	2709 (4.5)	36 (3.9)	17 (3.7)	ns	ns	ns
	Age older than 35 years	8354 (14.0)	186 (19.9)	99 (21.3)	< 0.0001	< 0.0001	ns
	Primiparity	22,153 (37.1)	455 (48.8)	208 (44.8)	< 0.0001	< 0.0001	ns
	Primigravidity	16,479 (27.6)	337 (36.1)	151 (32.5)	< 0.0001	0.02	0.03
	Single	21,456 (36.1)	352 (37.8)	165 (35.6)	ns	ns	ns
	High school or university	31,030 (54.2)	460 (51.5)	223 (53.1)	ns	ns	ns
	* Values are expressed as mean ± SD

Table 2

 2 Obstetric risk factors of the study population and comparison between EO-PE and LO-PE.

	Variables	No PE N = 59,665 N	LO-PE N = 933 N (%) EO-PE N = 464 N (%) LO-PE vs No PE p	EO-PE vs No PE p	EO-PE vs LO-PE p
		(%)			value	value	value
	Gestational diabetes	5978 (10.0)	132 (14.2)	55 (11.9)	0.004	NS	ns
	Pre-existing diabetes	816 (1.5)	29 (3.5)	18 (4.3)	< 0.0001	< 0.0001	ns
	Alcohol	389 (0.7)	6 (0.7)	3 (0.6)	ns	ns	ns
	Smoking during pregnancy	7397 (12.4)	83 (8.9)	42 (9.1)	0.002	0.04	ns
	Chronic hypertension	882 (1.5)	78 (8.4)	57 (12.3)	< 0.0001	< 0.0001	0.02
	Pre-pregnancy BMI (Kg/m 2 ) * 24.4 ± 5.8	26.9 ± 6.9	26.1 ± 6.0	< 0.0001	< 0.0001	ns
	BMI ≥ 30 kg/m 2	8975 (15.7)	252 (28.8)	99 (25.3)	< 0.0001	< 0.0001	ns
	IVF, ICSI	327 (0.5)	10 (1.1)	8 (1.7)	0.03	0.0007	ns
	Coagulopathy	171 (0.28)	5 (0.53)	6 (1.29)	ns	< 0.0001	ns
	Thyroid disease	783 (1.3)	15 (1.6)	12 (2.6)	ns	0.009	ns
	Renal disease	143 (0.24)	10 (1.07)	8 (1.7)	0.02	0.02	ns
	Hypercholesterolemia	26 (0.04)	3 (0.32)	2 (0.43)	0.01	0.01	ns

* Values are expressed as mean ± SD.

Table 3

 3 Delivery outcomes and complications of the study population and comparison between EO-PE and LO-PE.

	Variables	No PE N = 59,665 N	LO-PE N = 933	EO-PE N = 464	LO-PE vs No PE p	EO-PE vs No PE p	EO-PE vs LO-PE p	EO-PE vs LO-PE p
		(%)	N (%)	N (%)	value	value	value	adjusted for GA
	Caesarean section	9034 (15.1)	412 (44.2)	409 (87.5)	< 0.0001	< 0.0001	< 0.0001	< 0.0001
	Medically -Induced	12390 (20.7)	605 (64.8)	49 (10.6)	< 0.0001	< 0.0001	< 0.0001	< 0.0001
	delivery
	Abnormal Foetal Hearth	8731 (15.2)	244 (27.0)	157 (39.4)	< 0.0001	< 0.0001	< 0.0001	ns
	Rate
	Post-partum Haemorrhage 1406 (3.2)	40 (5.2)	5 (1.5)	0.002	ns	0.005	ns
	HELLP Syndrome	-	23 (2.5)	40 (8.6)	-	-	< 0.0001	ns
	Eclampsia	-	31 (3.3)	12 (2.6)	-	-	ns	ns

Table 4

 4 Foetal and neonatal outcomes in the study population and comparison between infants from EO-PE and LO-PE women.

	Variables	No PE N = 59,665 N	LO-PE N = 933 N	EO-PE N = 464 N	LO-PE vs No PE p	EO-PE vs No PE p	EO-PE vs LO-PE p	EO-PE vs LO-PE p
		(%)	(%)	(%)	value	value	value	adjusted for GA
	Gestational Age (weeks) *	38.4 ± 2.4	36.9 ± 1.8	30.4 ± 2.27	< 0.0001	< 0.0001	< 0.0001	na **
	Birth Weight (g) *	3085 ± 585	2617 ± 601	1228 ± 436	< 0.0001	< 0.0001	< 0.0001	na **
	Female sex	29,248 (49.0)	489 (52.4)	207 (54.5)	< 0.0001	< 0.0001	ns	ns
	SGA	4660 (7.8)	217 (23.3)	98 (21.1)	< 0.0001	< 0.0001	ns	ns
	Apgar@1 min ≤3	2709 (4.5)	93 (10.1)	111 (27.1)	< 0.0001	< 0.0001	< 0.0001	ns
	Need for intubation at	304 (1.8)	20 (2.8)	144 (35.1)	< 0.0001	< 0.0001	< 0.0001	0.007
	birth
	Intrauterine foetal death	524 (0.9)	12 (1.3)	31 (6.7)	ns	< 0.0001	< 0.0001	ns
	Medical termination of	279 (0.5)	1 (0.1)	20 (4.3)	ns	< 0.0001	< 0.0001	0.02
	pregnancy
	Early neonatal death	214 (0.4)	0 (0.0)	13 (2.8)	ns	< 0.0001	< 0.0001	ns
	Perinatal death	1017 (1.8)	13 (1.4)	64 (13.8)	ns	< 0.0001	< 0.0001	0.02