V S Afraimovich

V P Zhigulin

M I Rabinovich

On the origin of reproducible sequential activity in neural circuits

Robustness and reproducibility of sequential spatio-temporal responses is an essential feature of many neural circuits in sensory and motor systems of animals. The most common mathematical images of dynamical regimes in neural systems are fixed points, limit cycles, chaotic attractors, and continuous attractors (attractive manifolds of neutrally stable fixed points). These are not suitable for the description of reproducible transient sequential neural dynamics. In this paper we present the concept of a stable heteroclinic sequence (SHS), which is not an attractor. SHS opens the way for understanding and modeling of transient sequential activity in neural circuits. We show that this new mathematical object can be used to describe robust and reproducible sequential neural dynamics. Using the framework of a generalized high-dimensional Lotka-Volterra model, that describes the dynamics of firing rates in an inhibitory network, we present analytical results on the existence of the SHS in the phase space of the network. With the help of numerical simulations we confirm its robustness in presence of noise in spite of the transient nature of the corresponding trajectories. Finally, by referring to several recent neurobiological experiments, we discuss possible applications of this new concept to several problems in neuroscience.

Generation of sequential activity is one of the most important functions of neural circuits. For example, the high vocal center of songbirds produces a sequence of bursts that controls generation of unique and reproducible songs by the vocal system. The predator marine mollusk Clione produces a sequential hunting motion to find its prey and, in spite of the irregularity of this motion, the sequence of Clione's tail positions is also reproducible. Reproducible sequential neural activity is involved in sequential memory processes, fast recognition of stimuli by olfactory and auditory sensory systems and many other dynamical functions of the nervous system. What is the dynamical origin of such reproducibility? Which conditions do the synaptic connections of the neural circuits have to satisfy to generate reproducible transient sequences? In this paper we answer these questions using a generalized Lotka-Volterra model. We present a new mathematical construction-a stable heteroclinic sequence (SHS)-that describes reproducible sequence generation in the phase space of the dynamical model of a neural circuit. This construction keeps the trajectories in the vicinity of the heteroclinic skeleton due to the existence of a stable manifold on which the saddle points lie. These saddle points and the one-dimensional separatrices that connect them organize the SHS. We present rigorous results about the existence of the SHS that are supported by computer simulations of the dynamics along the SHS in presence of noise. We also discuss how the concept of SHS may help to understand observations of sequential dynamics in various neurobiological experiments.

I. INTRODUCTION

The principles of organization of neural circuits that are responsible for the formation of reproducible animal behavior are not completely known. The best understanding of such principles has been obtained for rhythmic motor functions like heart beating, breathing or walking. Such rhythmic functions are controlled by neural circuits that produce motor commands in the form of periodic spatio-temporal patterns of neural activity. The repetitive dynamics in these circuits is usually reproducible and can be described by a regular attractor, i.e., a limit cycle. However, many sensory, motor, and other neural systems demonstrate transient, nonperiodic activity, which, nevertheless, is reproducible from trial to trial. In these neural systems similar stimuli lead to the execution of similar sequential behaviors. The dynamics in neural circuits that control these finite-time nonperiodic behaviors is transient and, hence, cannot be implemented by a dynamical system with a limit cycle attractor.

In this paper we propose a new object in nonlinear dynamics of dissipative systems-a stable heteroclinic sequence (SHS), which represents a reproducible sequential pattern in the phase space of the system. We assume that the motor systems do not respond to all the detailed features of the output of sensory systems but just to the sequence of key events. Therefore, the reproducibility of spatio-temporal patterns that are generated by sensory neural circuits does not have to be perfect, but can be imperfect, i.e., the sequence of neuronal activations has to be reproducible, but their timing may not be. Such spatio-temporal sequential activity in neural ensembles is similar to the sequence of snapshots in a movie. It has been observed in the insect olfactory system, 1,2 mollusk sensory-motor system, 3 vocal tract of songbirds 4 and other neural systems. In those observations the time intervals between the "snapshots" were found to vary slightly from trial to trial but sequences of neuronal activations were found to be the same for the same stimulus and were changing with the stimulus. From the point of view of dynamical systems the observed reproducibility of sequential activity looks very intriguing because such transient dynamics is expected to depend strongly on initial conditions and fluctuations.

To explain this phenomenon we use as an example a Lotka-Volterra-type model of a neural network with inhibitory connections. Based on the results of computational experiments with this model we hypothesize that winnerless competition (WLC) 5 among the neurons in the network with dense inhibitory connections may explain the reproducibility of sequential neural activity. A WLC network of n neurons is a network which has a sequence of heteroclinic connections of saddle fixed points in its n-dimensional phase space. Each of these fixed points lies on a ͑n -1͒-dimensional stable manifold and these fixed points are connected along unstable directions by one-dimensional separatrices that induce transitions from one fixed point to the next. In this paper we formulate the necessary conditions on the connectivity of the WLC network that have to be satisfied in order for the network to exhibit reproducible sequential dynamics along the heteroclinic orbit.

Heteroclinic phenomena such as heteroclinic cycles and networks of several heteroclinic cycles have attracted the attention of many authors in the last 30 years (see, for example, Ref. 6). It has recently been proposed that the complex neural dynamics may be associated with these phenomena. 1 Although there are some results on the dynamics arising from the existence of heteroclinic contours and the behavior of the system in their vicinity, to the best of our knowledge results concerning dynamic behavior induced by an open heteroclinic sequence are absent. In our opinion, the SHS is the most suitable mathematical object for the description and analysis of reproducible sequential dynamics in biological and artificial neural networks and in complex systems in general.

The paper is organized as follows: In the second section we consider general features of the model, find conditions for the existence of a heteroclinic sequence, and formulate and prove a theorem about the reproducibility of sequential dynamics; in the third section we present results of numerical simulations that confirm our theoretical analysis and illustrate its generality; finally, we discuss the applicability of the results to several interesting problems of neuroscience.

II. ANALYTICAL RESULTS

A. The model

Our goal is to understand the reproducibility of sequential neural activity regardless of the exact timing of spikes. Hence, let us describe each neuron by one dynamical variable, a i ͑t͒, which corresponds to its average firing rate. Neurons in our model form a network with inhibitory synaptic connections. This model can also be extended to the case of excitatory-inhibitory networks, where excitatory principal neurons are connected through inhibitory interneurons like, for example, in the olfactory system of insects. 1 The dynamics of an inhibitory network of n = N + M neurons can be described by the following system of ordinary differential equations (ODEs):

da i ͑t͒ dt = a i ͑t͓͒ i ͑S͒ -͚ j=1 N+M ij ͑S͒a j ͑t͔͒ + i ͑t͒, ͑1͒
where N will be defined below [see inequality (32)], each a i ͑t͒ Ͼ 0 represents an instantaneous firing rate of ith neuron, i represents strength of stimulation of ith neuron, ij represents strength of inhibition of neuron i by neuron j and i ͑t͒ represents external noise. 5,7,8 In this model it is assumed that the stimulus S forms the matrix ij and increments i . Each increment i controls the time constant of an initial exponential growth from a i ͑t͒ =0.

Let us use a simplified model, 5 in which i Ͼ 0 when the stimulus S has a component at element i and i = 0 otherwise. If we single out only indices i for which i Ͼ 0 we obtain a system that could contain less than n equations, say N + M 1 , so (may be after relabeling of the indices) the system will have the form

da i ͑t͒ dt = a i ͑t͒ ͫ i -ͩ a i ͑t͒ + ͚ j i N+M 1 ij a j ͑t͒ ͪͬ , ͑2͒
where i Ͼ 0, i =1, ... ,N + M 1 and the external noise is assumed to be absent. Equation (2) implies, first of all, that the equilibrium at the origin O͑a i =0͒ is an unstable node.

B. Behavior of orbits around the origin

From now on we assume that

0 Ͻ i 1 Ͻ i , i = 1, ... ,N + M 1 , i i 1 . ͑3͒
It means that the strongly unstable manifold W O uu at O is the hyper plane a i 1 = 0 and the leading direction is the a i 1 -axis. In other words, all orbits go to O as t → -ϱ being tangent to the a i 1 -axis except for those belonging to W O uu . If one chooses randomly an initial point in a small neighborhood of O than almost surely the point will not belong to W O uu , and will leave the neighborhood following a i 1 direction.

C. The first saddle

Denote

by A 1 the equilibrium point ͑0, ... ,0, i 1 ,0, ... ,0͒ on the a i 1 -axis. By direct verification one can find the increments of the point A 1 ; they arei 1 and

͑ i -ii 1 i 1 ͒, i =1, ... ,N + M 1 , i i 1 . We assume that there ex- ists a number i 2 , 1ഛ i 2 ഛ N + M 1 such that i 2 -i 2 i 1 i 1 Ͼ 0 ͑4͒
and, in addition,

i -ii 1 i 1 Ͻ 0, 1 ഛ i ഛ N + M 1 , i ͕i 1 ,i 2 ͖. ͑5͒
It follows from these assumptions that the point A 1 is a saddle with one-dimensional unstable manifold. Moreover, the index i 2 indicates a direction in which the unstable separatrix of A 1 will go-see below.

For the sake of simplicity of the proof of further results, we assume that a i 1 -axis is a leading direction for the node A 1 on the hyperplane a i 2 = 0, i.e., we assume that

i -ii 1 i 1 Ͻ -i 1 , 1 ഛ i ഛ N + M 1 , i ͕i 1 ,i 2 ͖. ͑6͒
There is an important quantity, the so-called saddle value, that measures the evolution of a deviation from the stable and unstable manifolds along the orbits going through the neighborhood of A 1 . 9,10 A restriction of such a map to the plane P i 1 i 2 has the form

= c i 1 , ͑7͒
where c is a constant, is a deviation from the stable manifold, is a deviation from the unstable one and the saddle value is

i 1 = i 1 i 2 -i 2 i 1 i 1 . ͑8͒
If i 1 Ͼ 1 then the saddle A 1 is called dissipative: There is a contraction of a deviation after passing through a neighborhood A 1 .

D. The second saddle

The unstable direction of a point A 1 has a positive projection on the vector e i 2 = ͑0, ... ,0,1,0, ... ,0͒, where 1 corresponds to the i 2 th coordinate. Let us consider now the equilibrium point A 2 = ͑0, ... ,0, i 2 ,0, ... ,0͒. The increments at A 2 are -i 2 in the direction e i 2 and i -

ii 2 i 2 in the ith di- rection, i =1, ... ,N + M 2 , i i 2 .
We assume that there exists

i 3 ͕i 1 , i 2 ͖, 1ഛ i 3 ഛ N + M 1 such that i 3 -i 3 i 2 i 2 Ͼ 0 ͑9͒ and i -ii 2 i 2 Ͻ 0, ͑10͒ for 1 ഛ i ഛ N + M 1 , i ͕i 2 , i 3 ͖.
Under these assumptions the point A 2 is a saddle with one-dimensional unstable manifold. Moreover, the index i 3 shows the next point to be visited.

The intersection of hyperplanes

P i 1 i 2 = പ j i 1 ,i 2 N+M 1 ͕a j =0͖ is
the plane containing points A 1 , A 2 such that A 1 is a saddle point on P i 1 i 2 and A 2 is a stable node on it. The restriction of the system (2) on P i 1 i 2 has the form

a ˙i1 ͑t͒ = a i 1 ͑t͓͒ i 1 -͑a i 1 ͑t͒ + i 1 i 2 a i 2 ͑t͔͒͒, ͑11͒ a ˙i2 ͑t͒ = a i 2 ͑t͓͒ i 2 -͑a i 2 ͑t͒ + i 2 i 1 a i 1 ͑t͔͒͒, ͑12͒
and it follows from (4) and (10) that

i 2 -i 2 i 1 i 1 Ͼ 0 and i 1 -i 1 i 2 i 2 Ͻ 0. ͑13͒ This implies that if 1 -i 1 i 2 i 2 i 1 0 ͑14͒
then there are no equilibrium points in the region a i 1 Ͼ 0, a i 2 Ͼ 0. It is not difficult to show that the separatrix of the point A 1 belonging to the positive quadrant on P i 1 i 2 must go to attractor A 2 on the plane P i 1 i 2 , i.e., there is a heteroclinic connection ⌫ i 1 i 2 between points A 1 and A 2 . The increments of A 2 on P i 1 i 2 are -i 2 and ͑ i 1 -i 1 i 2 i 2 ͒, which is less than 0 because of (13). We assume that

-i 2 Ͻ i 1 -i 1 i 2 i 2 ͑15͒
or

i 1 i 2 Ͻ i 1 + i 2 i 2 .
Thus, the leading direction at A 2 is transversal to the a i 2 -axis on the plane

P i 1 i 2 .
It is natural to ask if this direction is the leading one on the hyperplane ͕a i 3 =0͖ (on which A 2 is a stable node). It will be so if all negative increments at A 2 will be less

than i 1 -i 1 i 2 i 2 , in other words if i -ii 2 i 2 Ͻ i 1 -i 1 i 2 i 2 ͑16͒
for i =1, ... ,N + M 1 , i ͕i 1 , i 2 , i 3 ͖. We assume that (15) and (16) are satisfied. A restriction of the local map in the neighborhood of A 2 on the plane P i 2 i 3 has the form

= c i 2 , ͑17͒
where c is a constant, is a deviation from the stable manifold, is a deviation from the unstable one and the saddle value is

i 2 = - i 1 -i 1 i 2 i 2 i 3 -i 3 i 2 i 2 = i 1 i 2 i 2 -i 1 i 3 -i 3 i 2 i 2 . ͑18͒
If i 2 Ͼ 1 then the saddle A 2 will be dissipative and an orbit comes closer to the unstable manifold of A 2 after going through a neighborhood of this point.

E. The mth saddle

Continuing inductively the procedure described above for the selection of saddle points and heteroclinic connections between them, we arrive at the mth level of the procedure. We suppose that the following assumptions were satisfied until now: (1) Selection of saddles.

The saddles A k = ͑0, ... ,0, i k ,0, ... ,0͒, k =2, ... ,m were selected in such a way that

i k+1 -i k+1 i k i k Ͼ 0, ͑19͒ i -ii k i k Ͻ 0, ͑20͒ for 1 ഛ i ഛ N + M 1 , i ͕i k+1 , i k ͖. Thus, all points A k ,2
ഛ k ഛ m are saddles with one-dimensional unstable manifolds (see Fig. 1). Moreover, the index i m+1 gives a way to proceed to the next step of induction. (2) Heteroclinic connections. Under the condition (which is supposed to be satisfied)

1 -i k-1 i k i k i k-1 0, ͑21͒
there is a heteroclinic orbit ⌫ i k-1 i k joining A k-1 and A k . This orbit belongs to the plane

P i k-1 i k = പ j i k-1 ,i k N+M 1 ͕a j =0͖.
This is easy to show taking into account the facts that there are no other equilibria on P i k-1 i k except for A k-1 and A k , the point A k is an attractor on P i k-1 i k , and there are no limit cycles on this plane. (3) Leading directions. Under the conditions (which are supposed to be satisfied)

-i k Ͻ i k-1 -i k-1 i k i k , ͑22͒ i -ii k i k Ͻ i k-1 -i k-1 i k i k , ͑23͒ for 1 ഛ i ഛ N + M 1 , i ͕i k-1 , i k , i k+1 ͖, the separatrix ⌫ i k-1 i k comes to
A k following a leading direction, transversal to the a i k -axis on the plane P i k-1 i k . Let us remark that the saddle value

i k = i k-1 i k i k -i k-1 i k+1 -i k+1 i k i k , ͑24͒
is defined for every saddle A k .

Consider now the point A m+1 = ͑0, ... ,0, i m+1 ,0, ... ,0͒. We assume that the following conditions are satisfied:

(1) A m+1 A k , 1ഛ k ഛ m. If A m+1 = A k then we have a heteroclinic contour. Such a situation has been studied in Ref. 11.

(2) There exists 1 ഛ i m+2 ഛ N + M 1 such that the following inequalities are satisfied

i m+2 -i m+2 i m+1 i m+1 Ͼ 0, ͑25͒ i -ii m+1 i m+1 Ͻ 0, ͑26͒ for 1 ഛ i ഛ N + M 1 , i ͕i m+1 , i m+2 ͖.
Thus A m+1 is a saddle with one-dimensional unstable manifold.

(3) Condition for heteroclinic connections is satisfied:

1 -i m i m+1 i m+1 i m 0. ͑27͒
It implies that there is a heteroclinic orbit ⌫ i m i m+1 joining A m and A m+1 . (4) Conditions for following the leading direction

-i m+1 Ͻ i m -i m i m+1 i m+1 , ͑28͒ i -ii m+1 i m+1 Ͻ i m -i m i m+1 i m+1 , ͑29͒ for 1 ഛ i ഛ N + M 1 , i ͕i m , i m+1 , i m+2 ͖ are supposed to be satisfied.
Now, under formulated conditions, there is a sequence of separatrices

⌫ m+1 ª ഫ k=1 m+1 A i k ഫ k=1 m ⌫ i k i k+1 , ͑30͒
joining A i 1 and A i m+1 . Let us assume that we can follow this inductive procedure until k = N, i.e., the conditions (19)-( 23) are satisfied for k =2, ... ,N. Thus, there exists a heteroclinic sequence ⌫ N consisting of saddle points A k and separatrices ⌫ i k i k+1 that join them. This sequence can be attracting in a sense. Let

1 ª i 1 , 2 ª i 1 • i 2 , ... , N ª ͟ k=1 N i k = i 1 i 2 -i 2 i 1 i 1 • ͟ k=2 N i k-1 i k i k -i k-1 i k+1 -i k+1 i k i k . ͑31͒
Then the following statement holds. Theorem: Assume that the conditions (19), (20), (22), and (23) are satisfied for k =1, ... ,N and the conditions (21) are satisfied for k =2, ... ,N. Assume also that k Ͼ 1, k = 1, ... ,N.

͑32͒

Then for any sufficiently small neighborhood U of the sequence ⌫ N there exists ⑀ 0 = ⑀ 0 ͑U͒ Ͼ 0 such that for every initial condition a 0 = ͑a 1 0 , ... ,a N+M 1 0 ͒ in a small neighborhood of A 1 with 0 Ͻ a i 2 0 = ⑀ ഛ ⑀ 0 , one has a͑t͒ U for 0 ഛ t ഛ T͑⑀͒, where a͑t͒ is a solution of (2) going through a 0 and T͑⑀͒ =-K log ⑀ + C where C is a constant independent of ⑀ and FIG. 1. A stable open heteroclinic sequence in a neural circuit with WLC. This is a mathematical image of reproducible nonrhythmic sequential activity. W i is a stable manifold of the ith saddle fixed point (heavy dots). The trajectories in the vicinity of the SHS represent sequences with different timings; the time intervals between switches are proportional to the square root of the distance from the SHS.

K = 1 i 2 -i 2 i 1 i 1 + ͚ k=2 N i k-1 i k+1 -i k+1 i k i k . ͑33͒ Moreover, a͑t͒ \ V͑ഫ k=1 N A i k ͒ (where V is the ⑀ 0 -neighborhood of equilibria A i k ) consists of N segments, and dist͑a͑t͒,⌫ i k i k+1 ͒ ഛ C 1 ⑀ k , ͑34͒
for values of t for which a͑t͒ belongs to the kth segment, where C 1 is a constant independent of ⑀.

Thus, an orbit starting in an ⑀-vicinity of the stable manifold of A 1 spends a long time in a neighborhood of the chain ⌫ N following heteroclinic connections ⌫ i k i k+1 . In other words, one may observe the same sequence of switching from one equilibrium to another for a large number of orbits. This qualitative statement is true even without conditions (22) and (23), but constants K, C, and C 1 might be different.

The proof of Theorem is obtained in the same way as the proof of Theorem 1 in Ref. 11.

Remark: One can see that the inequalities in the conditions of Theorem are the same as the ones in Theorem 1 in Ref. 11. The difference is that in Theorem 1 of Ref. 11 these conditions guarantee the stability of the heteroclinic contour as an attractor while here they only imply stability of the SHS during finite interval of time. Such kind of statements while being standard in, say, the theory of averaging is not so familiar in the dynamical systems theory. In other words, we construct the sequence of mappings of transversal to heteroclinic orbits sections along orbits close to these heteroclinic orbits absolutely in the same way as we did in Ref. 11. The conditions of Theorem imply that all these mappings are contractions. Thus, any orbit going through a point close to the first heteroclinic orbit will remain close to the SHS until it intersects the last transversal section to the unstable separatrix of the nth saddle point. In Theorem 1 of Ref. 11 the first and the last sections were the same-the heteroclinic orbits formed a closed heteroclinic contour-therefore, the corresponding Poincare mapping was a contraction. Here the first and the last sections are different, and one can guarantee the closeness of the SHS and an orbit only during a finite interval of time T͑⑀͒.

Corollary: Under the conditions of Theorem and provided that (3) is satisfied and given a neighborhood U as in Theorem, there exists an open set V + in the positive octant in a small neighborhood V of the origin such that every orbit going through V + comes to U at a point a 0 = ͑a 1 0 , ... ,a N+M 1 0 ͒ with 0 Ͻ a i 2 0 = ⑀ ഛ ⑀ 0 , and then follows the heteroclinic sequence ⌫ N staying in U for the time T͑⑀͒.

The proof of Corollary directly follows from the fact [assumption (3)] that most orbits leave V following the a i 1 -axis.

Remark: In reality V + can be very large, it can contain, say, 90% of all positive initial points in V. It explains reproducibility of sequences of switching ͕A i k ͖ while one repeats choosing points in the small neighborhood V of the origin and looks at the resulting orbits.

Summarizing the results of this section one can say that:

(1) Among saddle-type equilibria of the system (4) there are the ones with one-dimensional unstable separatrices and they can be found by using the inequalities (19) and (20). (2) Such saddle points can be connected by heteroclinic orbits that form ordered collections of saddles, so-called heteroclinic sequences. (3) A heteroclinic sequence can be stable and sufficient conditions for its stability are given in inequalities (32). ( 4) If these conditions are satisfied then an orbit through a point in a neighborhood of the first saddle of the heteroclinic sequence will stay in a vicinity of this sequence until its last saddle point. (5) The last statement holds for any choice of an initial point belonging to an open set in the phase space, hence reproducibility of spatio-temporal responses to similar stimuli is guaranteed.

III. NUMERICAL RESULTS

A. Simulation procedure

Let us represent conditions of Theorem in a more convenient form. First, the inequalities (20) and (22) are equivalent to

i k-1 i k Ͻ i k-1 i k Ͻ i k-1 i k + 1. ͑35͒
Furthermore, we assume that i k Ͼ 1, k =1, ... ,N. Then the conditions k Ͼ 1 will be automatically satisfied. Thus, the inequality

i k = i k-1 i k i k -i k-1 i k+1 -i k+1 i k i k Ͼ 1, ͑36͒
implies that

i k+1 i k Ͼ -i k-1 i k + i k-1 i k + i k+1 i k , ͑37͒
and, together with (35), it implies the inequality

i k+1 i k Ͼ i k+1 i k -1. ͑38͒
Combining the last result with (19) we obtain

i k+1 i k -1 Ͻ i k+1 i k Ͻ i k+1 i k . ͑39͒
Also, condition (23) may be rewritten as

ii k Ͼ i k-1 i k + i -i k-1 i k . ͑40͒
In order to check predictions of Corollary, we performed numerical simulations of the dynamics in a network of N = 50 neurons described by Eq. (2) with Gaussian white noise with mean= 0.02 and = 0.015 added to its right hand side. [START_REF]Integration was performed using the Runge-Kutta 6(5) scheme with variable time step. The absolute error was set to 10 -15 and the relative error was set to[END_REF] Increments i ͑i =1, ... ,N͒ were set to random numbers that were taken uniformly from the interval (5;10). Elements of the connectivity matrix ij were then chosen according to the following rules: Chaos, Vol. 14, No. 4, 2004 Sequential activity in neural circuits 1127

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 141.20.241.2 On: Tue, 14 Oct 2014 10:49:54

i k-1 i k = i k-1 i k + 0.5, ͑41͒ i k+1 i k = i k+1 i k -0.5, ͑42͒ ii k = i k-1 i k + i -i k-1 i k + 2, ͑43͒
that satisfy the conditions (35), (39), and (40).

The resulting system of equations was simulated 10 times, each time starting from a random initial condition within the hypercube R ͑0;0.2͒ 50 . In each simulation the sequence of crossings of the threshold a =4 by a i was recorded. Ten such networks with random sets of i were simulated and, as a result, 100 sequences were obtained, 10 for each network.

B. Results

Equations (41)-(43) satisfy the conditions of the Corollary. However, these conditions are not sufficient to observe reproducibility for arbitrary initial points. Hence, one cannot expect that all the networks with random i and strengths of connections which satisfy (41)-(43) will exhibit reproducible sequential dynamics. Nevertheless, most of the networks produced identical sequences in all 10 trials. In Fig. 2 we show an example of 10 color-coded time series for a network with a particular realization of i . Simulations of each trial were started from a different random initial condition. In this plot each neuron is represented by a different color and its level of activity is represented by the saturation of the color. This figure illustrates that the sequence of neuronal activations in the network is reproducible and the time intervals between switchings of neuronal activity are roughly the same.

In Fig. 3 we plot an example of the time series a i ͑t͒ during the initial stage of switching dynamics in the network. This plot illustrates the fact that the system "forgets" the initial conditions very quickly and from the very beginning starts to follow the heteroclinic sequence.

IV. DISCUSSION

The processing of sequential discrete events is important on every level of organization of neural systems: From one neuron to neural ensembles and microcircuits. One of the general dynamical principles that can be responsible for sequence generation in neural systems is the winnerless competition (WLC) among neurons. 1,5,8,[START_REF] Afraimovich | [END_REF] The mathematical image of the WLC dynamics is a trajectory in the vicinity of a stable heteroclinic sequence (SHS) in the phase space of the system. Such a sequence (see Fig. 1) is a chain of saddlenode fixed points connected by one-dimensional separatrices which retains nearby trajectories in its vicinity. In neural systems the WLC dynamics is realized by inhibitory connections between different neural clusters: Inhibitory neurons sequentially inhibit (turn off) or stop to inhibit (turn on) excitatory neurons. Flexibility of WLC dynamics may be provided by the dependence of the identity of participating neural clusters on the stimuli.

Sequence generation in the chain-or layer-like networks of neurons may result from feedforward wave-like propagation of spikes (like the waves in synfire chains [START_REF] Abeles | Corticonics[END_REF] ). In contrast, WLC dynamics does not need a specific spatial organization of the network. However, the image of a wave is a useful one, because in the case of WLC a "wave" of neural activity propagates in the network along the SHS. Such a wave is initiated by a stimulus. The speed of the sequential switching depends on the noise level. [START_REF] Afraimovich | [END_REF] The noise controls the distance between the realized trajectories of the system and the SHS. For trajectories that get closer to the SHS the time that the system spends near semi-stable states (saddle-nodes), i.e., the interval between switchings, becomes longer.

The mechanism of reproducibility of transient sequential neural activity that we have analyzed above is quite general and, as our preliminary results indicate, does not depend on the details of the neuronal model. We conjecture that saddle points in the phase space of the neural network can be replaced by saddle limit cycles [START_REF] Melbourne | [END_REF] or even chaotic sets that would describe neural activity in more detail, as in typical spiking or spiking-bursting models. This is extremely important for neurobiological applications because it may help to build a bridge between the concepts of neural competition and of neural synchronization. Formulation and proof of theorems like those presented in this paper for the more detailed neural models is a very interesting mathematical problem in itself. As our recent results indicate, 15 reproducible nonperiodic sequences can be generated even in neural networks with random inhibitory and excitatory connections given that the networks are far from the excitatory-inhibitory balance. Deep mathematical understanding of the spatiotemporal dynamics in random networks is also lacking.

Olfactory system may serve as one example of a neural system that generates transient, but reproducible from trial to trial sequences of neuronal activations. Complex intrinsic dynamics in the antennal lobe (AL) of insects transform static sensory stimuli into spatio-temporal patterns of neural activity. 1,2 Recently, new experimental results about the reproducibility of the transient spatio-temporal AL dynamics have been published. 2 In their experiments the authors presented bees with different odors and recorded neural activity in the ALs using calcium imaging. They have analyzed transient trajectories in the projection neuron activity space and found that the trajectories representing different trials of stimulation with the same odor were very similar. It was shown that after a time interval of about 800 ms different odors are represented in the phase space by different static attractors, i.e., after sufficiently long odor presentations the transient spatio-temporal patterns converge to different spatial patterns of activity. However, the authors emphasize that due to the reproducibility of the transient dynamics some odors were recognized in the early transient stage as soon as 300 ms after the onset of the odor presentation. It is highly likely that the transient trajectories observed in these experiments represent realizations of a SHS.

We would like to add here two more possible biological applications of the proposed mechanism of reproducible sequence generation. The first one is related to the observations of sequential dynamics exhibited by the place cells in hippocampus. 16 The dynamical principle for the sequential memory based on WLC has been suggested in Ref. 17. The essence of the idea is that the sequential memory is encoded in a multidimensional dynamical system with a heteroclinic sequence, where each saddle point represents an event of a sequence that is to be remembered. Once the state of the system approaches one fixed point representing a certain event, it is drawn along an unstable separatrix toward the next fixed point (event) and so on. Using the biologically motivated model of the sequential spatial memory it has been shown that, after learning, the system with WLC dynamics is capable of associative retrieval of pre-recorded sequence of patterns. 17 The second application is related to the dynamics of the high vocal center (HVC) in the brain of songbirds. It was observed that HVC dynamics is strictly sequential, with each projection neuron generating just one burst per sequence. 4 The sequence of bursts is highly reproducible because adult male birds produce precisely the same song in each performance. Taking into account the fact that the HVC neural circuit contains a lot of inhibitory interneurons, it is reasonable to hypothesize that its inhibitory connections are nonsymmetric and are organized in a way that is similar to the one described above.

FIG. 2 .

 2 FIG. 2. (Color). Time series of ten trials, simulations of each trial were started from a different random initial condition. Only one neuron is active at any given time, except for very short initial transient time (see Fig.3).

ACKNOWLEDGMENTS

This work was partially supported by NSF Grant No. EIA-0130708 and DOE Grant No. DE-FG03-96ER14592. V.S.A. was supported by CONACyT Grant No. 42765-F.


Camille Couffignal
email: camille.couffignal@inserm.fr

C Laouénan

C Burdet

A Foucrier

M Wolff

Armand - Lefevre

F Mentré

L Massias

Population pharmacokinetics of imipenem in critically ill patients with suspected ventilatorassociated pneumonia and evaluation of dosage regimens

Keywords: Imipenem, Population pharmacokinetics, critical care, Ventilator-associated pneumonia

Objectives: Significant alterations in the pharmacokinetics (PK) of antimicrobials have been reported in critically ill patients. We describe PK parameters of imipenem in intensive care unit (ICU) patients with suspected ventilator-associated pneumonia (VAP) and evaluate several dosage regimens.

Methods: This French multicentre, prospective, open-label study was conducted in ICU patients with a presumptive diagnosis of Gram-negative bacilli VAP who empirically received imipenem I.V. q8h.

Plasma imipenem concentrations were measured during the 4 th imipenem infusion using 6 samples (trough, 0.5, 1, 2, 5 and 8 hours). Data were analysed with a population approach using the SAEM algorithm in Monolix 4.2. A Monte Carlo simulation was performed to evaluate six dosage regimens: 500, 750 or 1000mg with administration q6h or q8h. The pharmacodynamic target was defined as the probability of achieving a fractional time (fT) above minimum inhibitory concentration (MIC) greater than 40%.

Results: Fifty-one patients were included in the PK analysis. Imipenem concentration data were best described by a two-compartment model with three covariates (creatinine clearance, total body weight and serum albumin). Estimated clearance (between-subject variability) was 13.2 L/h (38%), and estimated central volume 20.4 L (31%). At an MIC of 4 µg/mL, the probability of achieving 40% fT >MIC was 91.8%, for 0.5-h infusions of 750mg q6h, 86.0% of 1000mg q8h and 96.9% of 1000mg q6h.

Conclusion:

This population PK model accurately estimated imipenem concentrations in ICU patients. The simulation showed that for these patients the best dosage regimen of imipenem is 750mg q6h and not 1000mg q8h.

What is already known about this subject?

In critically ill patients, there are significant alterations in antimicrobials pharmacokinetics (PK) and high MICs bacteria. Pharmacodynamic target for imipenem, widely used in this specific population, is based on time over MIC, but few data is available for imipenem PK in these patients and dosage regimen have not been evaluated.

What this study adds?

Using a PK population approach, our study showed in critically ill patients a slight increased clearance and twice increased distribution volume of imipenem, compared to healthy patients. We also demonstrated that for 2-to-4 µg/mL MICs bacteria, a 750 mg q6h dosage regimen allowed to reach a 40% fractional time over MIC.

Introduction

Imipenem was the first licensed antibiotic of the carbapenem class and has been widely used for more than 30 years, for hospital-related infections caused by resistant Gram-negative bacilli. Due to its broad spectrum, imipenem is often prescribed for initial empirical treatment of ventilatorassociated pneumonia (VAP) in critically ill patients with risk factors for multidrug-resistant Gramnegative bacilli [START_REF]Guidelines for the management of adults with hospital-acquired, ventilator-associated, and healthcare-associated pneumonia[END_REF][START_REF] Joseph | The role of carbapenems in the treatment of severe nosocomial respiratory tract infections[END_REF]. It is a hydrophilic molecule characterised by a half-life (𝑡 1 2 ⁄ ) of one hour, low plasma protein binding (< 20%) and predominantly renal excretion unchanged close to 70% [START_REF] Balfour | Imipenem/cilastatin: an update of its antibacterial activity, pharmacokinetics and therapeutic efficacy in the treatment of serious infections[END_REF][START_REF] Wolff | Les carbapénèmes[END_REF]. In healthy subjects, the clearance is 12.1 L/h and the volume of distribution (Vd) of the central compartment is 9.6 L after 1000 mg every 8 hours (q8h) with 0.5-hour infusion [START_REF] Lee | Comparison of 30-min and 3-h infusion regimens for imipenem/cilastatin and for meropenem evaluated by Monte Carlo simulation[END_REF]. Imipenem has a time-dependent bactericidal activity and the pharmacodynamic (PD) parameter associated with its bactericidal effect is the fractional time (fT) when concentration is above the minimum inhibitory concentration (MIC). Imipenem has a post-antibiotic effect of 2 to 6 hours against most Gramnegative bacilli. Antimicrobial activity is optimised when the fraction of time above MIC (fT > MIC) is greater than 40% [START_REF] Rodloff | Two decades of imipenem therapy[END_REF][START_REF] Mouton | Comparative pharmacokinetics of the carbapenems: clinical implications[END_REF], but for critically ill patients some studies suggest an optimal fT > MIC of 100% [START_REF] Craig | Pharmacokinetic/pharmacodynamic parameters: rationale for antibacterial dosing of mice and men[END_REF][START_REF] Zelenitsky | Pharmacodynamics of empirical antibiotic monotherapies for an intensive care unit (ICU) population based on Canadian surveillance data[END_REF].

In critically ill patients , the pharmacokinetic (PK) properties of drugs are modified with an increase in Vd, fluctuation of plasma clearance, presence of oedema, and drug-drug interaction [START_REF] Roberts | Pharmacokinetic issues for antibiotics in the critically ill patient[END_REF][START_REF] Pea | Antimicrobial therapy in critically ill patients: a review of pathophysiological conditions responsible for altered disposition and pharmacokinetic variability[END_REF][START_REF] Mehrotra | Antibiotic pharmacokinetic and pharmacodynamic considerations in critical illness[END_REF][START_REF] Boucher | Pharmacokinetic changes in critical illness[END_REF], resulting in a lesser or higher drug exposure.

In addition to changes in PK characteristics of ICU patients, there is a worrisome increase in the incidence multiresistant Gram-negative bacilli, especially in the ICU. In this context, dosage regimens of antibiotics in the ICU must be adapted. Currently clinicians tend to increase the doses of antibiotics or change the dosage schedule without customizing antibiotic regimens according to the host and the offending pathogen.

Although imipenem is widely used in critically ill patients, data allowing optimisation of its administration are surprisingly scarce. Published PK studies with data in these patients are either of imperfect design [START_REF] Belzberg | Imipenem levels are not predictable in the critically ill patient[END_REF] or have a small number of subjects [START_REF] Mckindley | Pharmacokinetics of aztreonam and imipenem in critically ill patients with pneumonia[END_REF][START_REF] Tegeder | Tissue distribution of imipenem in critically ill patients[END_REF][START_REF] Novelli | Pharmacokinetic evaluation of meropenem and imipenem in critically ill patients with sepsis[END_REF][START_REF] Sakka | Population pharmacokinetics and pharmacodynamics of continuous versus short-term infusion of imipenem-cilastatin in critically ill patients in a randomized, controlled trial[END_REF][START_REF] Jaruratanasirikul | Comparison of the Pharmacodynamics of Imipenem in Patients with Ventilator-Associated Pneumonia Following Administration by 2 or 0.5 H Infusion[END_REF]. Among these studies, three have specifically analysed imipenem PK variability in ICU patients with VAP.

Some authors have evaluated several dosage regimens of antimicrobials in critically ill patients. For aminoglycoside antibiotics, Conil et al. [START_REF] Conil | Tobramycin disposition in ICU patients receiving a once daily regimen: population approach and dosage simulations[END_REF] showed the impact of adapting the regimen on PD targets (80 < AUC < 125 mg.L -1 .h and peak > 10 mg.L -1 ) after simulation dosage regimens in ICU patients with nosocomial infections. For meropenem, an antibiotic of the same class as imipenem, Crandon et al [START_REF] Crandon | Optimization of meropenem dosage in the critically ill population based on renal function[END_REF] evaluated the concentration-time profile in ICU patients with VAP in order to limit the potential inadequacies noted for current dosage regimens.

The aim of our study was to estimate the PK parameters of imipenem and their variability in ICU patients with suspected VAP, using a population approach to determine the influence of clinical and biological covariates for imipenem. We performed a Monte Carlo simulation to evaluate several dosage regimens based on the PD parameter (fT > MIC) for the range of clinical relevant MICs in ICU.

Materials and Methods

Study design and population

Sampling Procedure and Analytical Methods

Imipenem concentrations were measured at steady state after the 4 th dose i.e. between 24 and 32 hours after the first infusion of imipenem. Six blood samples per patient were collected immediately before and at 0.5, 1, 2, 5 and 8 hours after the 4 th infusion for concentration measurement.

Blood samples were retrieved from 4 mL of heparin and immediately centrifuged at 4000 rpm.

Plasma was then stabilised within ½ hour after collection, by 4-morpholine propane sulphonic acid (MOPS) in ethylene glycol and immediately frozen at -80 °C. Plasma imipenem concentrations were determined after processing the samples by ultrafiltration, using high-performance liquid chromatography (HPLC) on an Interchrome © YP5C18 25QS reverse phase column (length 25 cm, internal diameter 4.6 cm). UV detection was performed at 302 nm [START_REF] Garcia-Capdevila | Determination of imipenem in plasma by high-performance liquid chromatography for pharmacokinetic studies in patients[END_REF]. Chromatographic peaks were integrated and imipenem concentrations calculated using Empower 2 software Water®. The lower limit of quantification (LOQ) was 0.5 mg/L.

Blood sample analysis was centralised in the pharmacology-toxicology laboratory of the Hôpital Bichat, AP-HP, Paris, France.

Population pharmacokinetic model building

Population PK analysis was performed using MONOLIX 4.1.2 software (www.lixoft.eu). Population PK parameters were estimated by maximum likelihood using the stochastic approximation expectation maximisation (SAEM) algorithm [START_REF] Lavielle | Estimation of Population Pharmacokinetic Parameters of Saquinavir in HIV Patients with the MONOLIX Software[END_REF]. The SAEM algorithm is an expectation maximisation (EM) algorithm extension in the nonlinear mixed-effects models where the parameter estimation was computed by the maximum likelihood estimator of the parameters without any approximation of the model as linearization. Briefly, SAEM converges to maximum likelihood estimates by repeatedly alternating between the E and M steps. Then, the expectation of the complete likelihood is computed according to a stochastic approximation [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF].

The full maximum likelihood estimation allows to take into account the data below quantification limit (BQL) [START_REF] Samson | Extension of the SAEM algorithm to left-censored data in nonlinear mixed-effects model: Application to HIV dynamics model[END_REF]. BQL data are considered as left-censored observations, indeed in that case the data 𝑦 𝑖𝑗 is not observed but we only know that it is below the LOQ. The extension of the SAEM algorithm in MONOLIX to consider BQL realized a simulation of the left-censored data in a right-truncated Gaussian distribution with an integration below limit of quantification to obtain probability of BQL. It is very similar to the method call 'M3' in NONMEM for handling BQL data [START_REF] Beal | Ways to fit a PK model with some data below the quantification limit[END_REF].

Structural and statistical model

In the first step, a basic population PK model without covariates was developed. For the structural PK model, one-and two-compartment models were compared. Exponential random effects were assumed to describe between-subject variability: e.g. for clearance (CL) of subject i, 𝐶𝐿 𝑖 = 𝐶𝐿 𝑝𝑜𝑝 × e 𝜂 𝐶𝐿,𝑖 where 𝐶𝐿 𝑝𝑜𝑝 is the population parameter estimate, 𝜂 𝐶𝐿,𝑖 is the individual random effect. The random effects were first supposed to be independent with diagonal variance-covariance matrix Ω and then possible correlations between random effects were tested in this variance-covariance matrix. Additive, proportional and combined error models were tested. The most appropriate pharmacostatistical model was selected on the basis of the following criteria: (i) smaller value of Bayesian information criterion (BIC); (ii) adequate goodness-of-fit (GOF) plots; (iii) low relative standard error (RSE) in estimated PK parameters.

Covariate analysis

From the basic model, twelve covariates were studied and chosen for their impact on the PK parameters specifically in the ICU in accordance with published data. These 12 covariates were: age, gender, total body weight at inclusion and total body weight change (between the 4 th dose and admission); three specific ICU scores, namely SAPS II [START_REF] Gall | A new Simplified Acute Physiology Score (SAPS II) based on a European/North American multicenter study[END_REF], the SOFA score [START_REF] Vincent | The SOFA (Sepsisrelated Organ Failure Assessment) score to describe organ dysfunction/failure. On behalf of the Working Group on Sepsis-Related Problems of the European Society of Intensive Care Medicine[END_REF] and the oedema score (ES) [START_REF] Gómez | Altered pharmacokinetics of ceftazidime in critically ill patients[END_REF]; serum albumin and four-hour creatinine clearance (CrCl4h) [START_REF] Pickering | Four hour creatinine clearance is better than plasma creatinine for monitoring renal function in critically ill patients[END_REF]; positive end-expiratory pressure (PEEP), PaO2/FiO2 ratio, and the presence of septic shock. These covariates were recorded at the 4 th dose of imipenem except for SAPS II and weight, which were measured both at admission and inclusion. Urine samples for CrCl4h were collected when the 4 th infusion of imipenem had started.

CrCl measurement over 4 hours was assumed to be a true reflection of renal function during the 4 th infusion [START_REF] Herrera-Gutiérrez | Replacement of 24-h creatinine clearance by 2-h creatinine clearance in intensive care unit patients: a single-center study[END_REF][START_REF] Herget-Rosenthal | How does late nephrological co-management impact chronic kidney disease? -An observational study[END_REF]. Missing values for tested covariates were imputed to the median value observed in the analysis population.

The parameter-covariate relationships were modelled multiplicatively as follows (e.g. for imipenem clearance CL): for continuous covariates, 𝐶𝐿 𝑖 = 𝐶𝐿 𝑝𝑜𝑝 × ( 𝐶𝑂𝑉 𝑖 𝐶𝑂𝑉 𝑚𝑒𝑑𝑖𝑎𝑛 ) 𝛽 × e 𝜂 𝐶𝐿,𝑖 where 𝛽 is the covariate effect to be estimated, 𝐶𝑂𝑉 𝑖 is the value for the subject 𝑖; 𝐶𝑂𝑉 𝑚𝑒𝑑𝑖𝑎𝑛 is the median value of covariates; for binary covariates, 𝐶𝐿 𝑖 = 𝐶𝐿 𝑝𝑜𝑝 × 𝑒 𝛽.𝐶𝑂𝑉 𝑖 × e 𝜂 𝐶𝐿,𝑖 where 𝐶𝑂𝑉 𝑖 takes a value of 0 or 1. For all covariates, binary or continuous, the unit of 𝛽 is the log of the unit of the associated parameter.

Covariates were selected with a forward method using BIC [START_REF] Mould | Basic concepts in population modeling, simulation, and model-based drug development-part 2: introduction to pharmacokinetic modeling methods[END_REF]. First, a model with one covariate was selected with the smallest BIC. Then, the model with two covariates was selected similarly. The addition of covariates was stopped when no further decrease of BIC was obtained. The covariates model was finalised with a backward selection, removing covariates one by one, using the Likelihood Ratio Test (LRT). A covariate was retained in the model if the LRT was significant (p < 0.05) when it was removed from the full model. In the final model, the 95% confidence interval of each parameters was determined from 1000 nonparametric bootstraps based resampling [START_REF] Thai | Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed-effects models: a simulation study in population pharmacokinetics[END_REF].

Model evaluation

Evaluation of the model was based on GOF plots. The model was first evaluated using observations versus individual and population predictions plots and usual residual-based plots (individual weighted residuals [IWRES] plot and population weighted residuals [PWRES] plot). It was then assessed using simulation-based plots, (visual predictive check (VPC) plot and normalised prediction distribution error (NPDE) versus time. The VPC plot showed the 10 th , 50 th , and 90 th percentiles of observed data over time and their corresponding 90% prediction intervals calculated from 500 Monte Carlo samples (simulated using the model, the parameter estimates and the design of the dataset). NPDE was built from the percentile derived from VPC prediction. The plot of NPDE takes into account the full predictive distribution of each individual observation and the various imipenem doses. As only few patients had different doses of imipenem and as we plot NPDE, we did not perform a prediction corrected VPC [START_REF] Bergstrand | Prediction-corrected visual predictive checks for diagnosing nonlinear mixed-effects models[END_REF].

Model evaluation was performed for both the basic model and the final model with covariates.

Monte Carlo simulation for dosage regimen evaluation

A Monte Carlo simulation was performed using the final PK model with covariates to predict the distribution of plasma imipenem concentrations and to estimate the PD parameter fT > MIC for several current dosage regimens and various MIC values. Six usual dosage regimens were studied: 500, 750 and 1000 mg with administration q6h or q8h. We simulated 1000 patients with a set of covariates re-sampled among the observed covariates of included patients and a vector of random effects drawn from the estimated distribution. The concentration-time profile of the 1000 virtual patients was simulated at steady state for the six dosage regimens.

The MIC targets were selected from the European Committee on Antimicrobial Susceptibility Testing (EUCAST [START_REF] Eucast | Breakpoint table for interpretation of MICs and zone diameters[END_REF]) data and ranged from 0.06 to 32 µg/mL. Two specific MIC, 2 µg/mL and 4 µg/mL, were studied. These MICs were the limited sensitivity breakpoint of imipenem currently observed for Gram-negative bacilli isolated in the ICU (Enterobacteriaceae species and Pseudomonas aeruginosa, respectively).

The time for which the imipenem concentration remained above the MIC at steady state was calculated as a cumulative percentage over a 24-hour period and the probability of pharmacodynamic target attainment (PTA) was assessed as a fraction that achieved 40% fT > MIC or 100% fT > MIC.

Results

Patients

Sixty-three patients were included in the IMPACT study. Twelve patients were excluded from the PK analysis: three lacking a kinetic profile and nine who did not receive four doses of imipenem. Fiftyone patients were included in the PK analysis, 41 of whom were males (80%), ranging in age from 28 to 84 years (median 60 years). At inclusion, median total body weight was 77 kg (range [45-126]). All patient characteristics are summarised in Table 1. Reasons for admission to the ICU were medical in 40 patients (78%) and surgical for 11 patients (22%), and the SAPS II at admission was 43 . The median duration of stay in the ICU and of mechanical ventilation before inclusion was 8 days [1-60] and 8 days , respectively. Antibiotic therapy was prescribed to 48 patients (94%) in the three months before admission, including 11 patients (30%) who previously received imipenem.

Four patients (9%) received 500 mg of imipenem, 15 (29%) 750 mg and 32 (62%) 1000 mg with the same dose interval q8h.

Population pharmacokinetic analysis

A total of 297 samples were available for PK modelling with a median of 6 samples [START_REF] Balfour | Imipenem/cilastatin: an update of its antibacterial activity, pharmacokinetics and therapeutic efficacy in the treatment of serious infections[END_REF][START_REF] Wolff | Les carbapénèmes[END_REF][START_REF] Lee | Comparison of 30-min and 3-h infusion regimens for imipenem/cilastatin and for meropenem evaluated by Monte Carlo simulation[END_REF][START_REF] Rodloff | Two decades of imipenem therapy[END_REF] per individual (Figure 1). Imipenem concentrations at peak (0.5 h) and trough were 34.1 mg/L [12.3-67.5] and 1.9 mg/L [0.5-10.1], respectively. Nine percent of imipenem concentrations were below the limit of quantification (BQL). One patient received the 4 th dose 5 hours late.

Imipenem PK concentrations were best described by a two-compartment model. An exponential random effects model described the between-subject variability in clearance CL and volume of distribution of the central compartment V1. Since the variability of intercompartmental clearance Q and the volume of distribution of the peripheral compartment V2 were very low, the betweensubject variability was not estimated and was taken as zero. A proportional model was used to describe the residual variability.

As shown in Table 2, estimated imipenem CL was 13 L/h, Q 10.1 L/h and the volumes V1 and V2 were 22.4 L and 9.9 L, respectively. A correlation between CL and V1 was retained in the basic model and estimated as 0.48. The GOF plots of the basic model were satisfactory (plots not shown).

Model with covariates

The best model with one covariate included the effect of 4-hour creatinine clearance (CrCl4h) on CL.

Covariate selection was continued up to a model with four covariates; the model with five covariates had a larger BIC (Table 3).

The backward selection was then performed from the model with the four following covariates: CrCl4h, age, serum albumin and total body weight. Only three covariates were significant using the LRT and kept in the final model: CrCl4h on CL, serum albumin (imputed to median value for 8 patients with missing data) and total body weight on V1. Imipenem CL was found to increase with CrCl4h. V1 was found to increase with total body weight and decrease with serum albumin (Figure 2).

The introduction of CrCl4h alone reduced the variability of CL (ωCL) from 48% to 38%. The introduction of weight and serum albumin reduced the variability of V1 from 48% to 31%. The final PK parameters are summarised in Table 2. All were reliably estimated, as reflected by the small RSEs from observed Fisher Information Matrix. The results of bootstrap medians and 95% CI were consistent except for the between-subject variability ωV1 and the correlation. Nevertheless the bootstrap analysis confirmed the reliability and robustness of the parameter estimates and thus the final model with covariates was representive. Estimated parameters were similar in the analysis of the 43 patients with no missing albumin data (results not shown).

Model evaluation

The GOF plots of the final PK model with covariates are shown in Figure 3 and

Monte Carlo simulation for dosage regimen evaluation

Using the simulated concentration-time profiles at steady-state, the PTA (40% or 100% fT > MIC) was calculated for the current dosage regimens 500, 750 and 1000 mg q6h or q8h. As shown in Figure 5 (a), all simulated patients had a fT > MIC greater than 40% for MIC from 0.06 to 1 µg/mL for the 6 dosage regimens. For MIC = 2 µg/mL, 86% of patients had the PTA at 40% with 500 mg q8h, 96.9% with 500 mg q6h; 95.3% with 750 mg q8h, 99.1% with 750 mg q6h; 97.9% with 1000 mg q8h and 99.4% with 1000 mg q6h. Figure 5 (b) shows the probability of fT > MIC greater than 100% with the 6 different dosage regimens. The percentage of patients was higher with the q6h regimen than with the q8h regimen, whatever the dose. For MIC = 4 µg/mL, 5% of simulated patients had the PTA at 100% with 500 mg q8h and 18.7% with 500 mg q6h; 14.3% with 750 mg q8h and 32.5% with 750 mg q6h; 20.9% with 1000 mg q8h and 45% with 1000 mg q6h. These results were confirmed by the simulated median concentration-time profile after four doses of imipenem (study protocol) as shown in Figure 6. The median patient with the 1000 mg q8h regimen did not achieve the PD target of MIC = 2 µg/mL. The median patient with 750 mg or 1000 mg q6h achieved the PD target of MIC = 2 µg/mL, but not for MIC = 4 µg/mL. We also explored the impact of each of the three significant covariates, namely CrCl4h, total body weight and serum albumin, on PTA. In Table S1 (of supplementary materials), we computed the PTA for the dose of 1000 mg q8h (the dosage regimen of the protocol) and 750 mg q6h (same daily dose) for three percentiles of each covariate (10 th , 50 th and 90 th ), assuming the two remaining covariates were at their median value. For both dosing regimens, 40% fT > MIC was obtained for the target MICs 2 µg/mL and 4 µg/mL for all values of the covariates. Figure S1 shows the concentration profile at steady state for 750 mg q6h for the various covariate values. It illustrates the rather limited impact of covariates on fT > MIC for that dosage regimen. median CrCL = 54.1 mL/min, and were lower than those estimated in our population. In our study, we observed high values of creatinine clearance. These high values are consistent with the hyperdynamic state of sepsis patients and confirmed the physiological impact on the PK parameters [START_REF] De Paepe | Pharmacokinetic and pharmacodynamic considerations when treating patients with sepsis and septic shock[END_REF]. It is therefore necessary to regularly control this creatinine clearance parameter in ICU to limit an effect on the clearance of imipenem although this effect was limited in our final PK model (Figure S1). A similar process of recorded creatinine clearance was previously described by Belzberg et al. ( 14) during 2 hours with maximum value of 408 mL/min. in accordance with our study, this study did not found an influence of creatinine clearance on imipenem clearance but an increase of distribution volume.

The hydrophilic nature of imipenem makes it sensitive to changes in the distribution of body fluids.

Its volume of distribution is affected by all disorders resulting in an increase of the extracellular compartment such as sepsis or clinically revealed by edema. In our study, we did not find any effect of edema on imipenem volume of distribution in the PK model but total body weight and serum albumin were found to influence significantly the distribution of imipenem and probably were reflected the physiological characteristics encountered in the ICU. The same increase of central volume was also observed for another antibiotic class, the aminoglycosides. Tanigawara et al. [START_REF] Tanigawara | Population pharmacokinetics of Arbekacin in patients infected with methicillin-resistant Staphylococcus aureus[END_REF] showed a significant increase of volume of distribution in a comparative study between healthy subjects and patients with pneumonia or sepsis treated by arbekacin hydrophilic and low protein binding antibacterial agents as imipenem. For results of Figure 2 and total body weight estimate coefficient, we evaluated the volume with a coefficient of total body weight fixed to one. The volumes expressed as L/kg were very similar between the final model and the total body weight coefficient model, 0.26 L/kg and 0.27 L/kg respectively. Due to the hydrophilic property of this antibacterial agent, we wished to evaluate other weight metrics (as ideal body weight or lean body weight) but these parameters could unfortunately not be collected during patient monitoring. No other covariates, especially ICU scores, were found in the PK model and considered for determination of the dosage regimen of imipenem.

With the PK results, we also performed a Monte Carlo simulation to evaluate several dosage regimens with doses given q6h or q8h. In the context of suspected VAP due to Gram-negative bacilli, we focused our simulations on the target MICs of 2 and 4 µg/mL (sensitivity breakpoint of Enterobacteriaceae species and Pseudomonas aeruginosa, respectively). With the same daily dose of 3 g, a q6h infusion led to a PD objective greater than 40% fT > MIC for those two target MICs. Despite their impact on the variability of PK parameters, the covariates lead to rather small changes in PK parameters and concentration profiles and thus have a limited effect on PTA, and we show rather small changes in the PTA 40% fT > MIC. Our evaluation of the 3 g daily dose has confirmed that q6h is a good dosage regimen for use in the ICU. This dose did not exceed the threshold of toxicity and the q6h regimen was optimised to take into account the higher PK variability seen in critically ill patients.

Conclusion

Our results demonstrate that imipenem pharmacokinetics vary in ICU patients. Imipenem clearance CL and central volume V1 were best estimated with three covariates whose influence on pharmacokinetic estimates was limited. Using population pharmacokinetic parameters, we showed that an infusion 750 mg q6h dosage regimen (3 g daily dose) is needed to achieve adequate pharmacodynamics, i.e., a fraction of time above MIC greater than 40% for usual the MICs of 2 and 4 µg/mL. ) -1.1 .

Figure legends:

Figure 4 .

 4 The model adequately described the observations as shown by the plots of observations versus population and individual predictions with the exception of the highest concentrations. Moreover, the NPDE plot versus predictions and the IWRES show no trend. The VPC plot and the NPDE plot presented in Figure 4 as a function of time from first dose indicate a good predictive performance of the model.

Figure 1 :

 1 Figure 1: Spaghetti plot of imipenem concentrations versus time following four doses for the 51 ICU patients included in the analysis. Data above LOQ are presented as blues circles, BQL data as red circles at LOQ.

Figure 2 :

 2 Figure 2: Relationship between estimated individual pharmacokinetic parameters and covariates: (a) clearance vs. 4-hour creatinine clearance (b) central volume vs. total body weight (c) central volume vs. serum albumin. Model predictions are displayed as the red curve.

Figure 3 :

 3 Figure 3: Goodness-of-fit plots for the final model with covariates: (a) Observations and (b) population weighted residuals (PWRES) versus population predicted values; (c) observations and (d) individual weighted residuals (IWRES) versus individual predicted values. Observations are plotted as blue circles and BQL data as red circles. LOWESS smoothed curve are plotted as blue curves.

Figure 4 :

 4 Figure 4: (a) Visual Predictive Check (VPC) and (b) Normalized Prediction Distribution Error (NPDE) versus time since first dose for the final model. VPC details: the solid green lines indicate the 10 th, 50 th and 90 th percentiles for observed data. The shaded blue and pink areas represent 90% prediction intervals from the corresponding percentiles calculated from simulated data. Observations are plotted as blue circles and BQL data as red circles.

Figure 5 :

 5 Figure 5: Simulated probabilities of pharmacodynamic target attainment versus MIC for various imipenem current dosage regimens at steady state (a) 40 % fT > MIC and (b) 100 % fT > MIC. Vertical lines are displayed for MIC = 2 µg/mL and MIC = 4 µg/mL which are the thresholds currently observed for Gram-negative bacteria in the ICU.

Figure 6 :

 6 Figure 6: Predicted concentrations of imipenem for median value of parameters for 1000 mg q8h or 1000 mg q6h or 750 mg q6h.

*

  10 th = 17 mL/min; 50 th = 86.4 mL/min and 90 th = 258 mL/min ** 10 th = 53 kg; 50 th = 77 kg and 90 th = 111 kg *** 10 th = 11 g/L; 50 th = 18 g/L and 90 th = 23 g/L

Figure S1 :

 S1 Figure S1 : Predicted steady-state concentrations of imipenem for 750 mg q6h dosage regimen with percentile values (10th, 50th and 90th) of the three significant covariates: With (a) Creatinine clearance percentiles: 10th = 17 mL/min; 50th = 86.4 mL/min; 90th = 258 mL/min (b) Body weight: 10th = 53 kg; 50th = 77 kg; 90th = 111 kg (c) Serum albumin: 10th = 11 g/L; 50th = 18 g/L; 90th = 23 g/L Vertical lines are displayed for MIC = 2 μg/mL and MIC = 4 μg/mL.

Table 1 : Characteristics at inclusion or at time of PK collection (4 th dose) of the 51 ICU patients included in the PK analysis

 1

	Parameters	Value*
	At inclusion
	Male	41 (80 %)
	Age (years)	60 [28-84]
	Total body weight (kg)	77 [45-126]
	SAPS II	40 [19-74]
	At time of 4 th dose
	Weight change # (kg)	1.1 [-18.1-19.1]
	SOFA	6 [2-14]
	Oedema score	7 [0-18]
	Serum albumin (g/L)**	18 [10-28]
	CrCl4h (mL/min)	86.4 [9.1-571.4]
	Shock	18 (35 %)
	PEEP (cmH2O)	6 [0-13]
	PaO2/FiO2	182 [81-346]

* Values are expressed as median [min-max] or number (percent)

# between the 4 th dose and admission ** median value for 9 patients

Table 2 : Population PK parameters of imipenem in 51 ICU patients

 2

		Basic model			Final model
		Value	RSE * (%)	Value	RSE * (%)	p-value **	Median bootstrap #	95% CIs bootstrap #
	Fixed effects
	CL (L/h)	13.0	6	13.2	5			13.2	11.4 -15.3
	βCrCL4h (log L/h)	-	-	0.2	19	6.4 10 -5	0.25	0.1 -0.4
	V1 (L)	22.4	9	20.4	7			19.8	14.9 -25.4
	βWeight (log L)	-	-	1.3	17	1.3 10 -4	1.2	0.6 -2.2
	βSerum albumin (log L)	-	-	-1.1	18	1.8 10 -4	-1.0	-1.8 --0.5
	Q (L/h)	10.1	28	12.2	25			12.3	4.7 -20.3
	V2 (L)	9.9	14	9.8	13			10.5	6.9 -13.7
	Between-subject variability
	ωCL (%)	48	10	38	13			36	26 -49
	ωV1 (%)	48	15	31	18			22	1 -45
	Correlation
	𝜂 𝐶𝐿 𝑖 𝜂 𝑉 1 𝑖	0.48	29	0.51	28			0.79	-1 -1
	Residual variability
	σ (%)	33	4	33	3			34	26 -41
	BIC	1595	-	1560	-			-
	𝑎𝑙𝑏𝑢𝑚𝑖𝑛 ( 18					𝐶𝑟𝐶𝐿 𝑖 86.4	) 0.2 𝑎𝑛𝑑 𝑉 1𝑖 = 20.4 × ( 𝑤𝑒𝑖𝑔ℎ𝑡 77	) 1.3 ×

* RSE: relative standard error; ** Likelihood ratio test (LRT); # from 1000 bootstrap resampling.

Final population PK covariate model is: 𝐶𝐿 𝑖 = 13.2 × (

Table S1 : Expected fractional time above MIC (fT > MIC) for two target MICs 2 and 4 µg/mL for 1000 mg q8h and 750 mg q6h dosage regimen and population parameters for the 10 th , 50 th and 90 th percentiles of the three significant covariates.

 S1

			MIC = 2 µg/mL			MIC = 4 µg/mL
	1000 mg q8h	10 th	50 th	90 th	10 th	50 th	90 th
	CrCl4h*	100	86.5	66.5	97.2	63.7	48.1
	Body weight**	72.7	86.5	100.0	53	63.7	77.4
	Serum albumin***	100.0	86.5	78.0	80.0	63.7	57.7
			MIC = 2 µg/mL			MIC = 4 µg/mL
	750 mg q6h	10 th	50 th	90 th	10 th	50 th	90 th
	CrCl4h*	100	100	79.9	100	74.8	55.2
	Body weight**	88.0	100	100	63.1	74.9	92.0
	Serum albumin***	100	100	94.7	95.3	74.9	67.7

Funding Statement

This study was funded by the Contrat d'Initiation à la Recherche Clinique 2006 (Assistance Publique-Hôpitaux de Paris, Département de la Recherche Clinique et du Développement, CRC 06049)

Discussion

We studied the pharmacokinetics of imipenem after I.V. infusion in 51 critically ill patients hospitalised in an ICU with suspected Gram negative VAP, using a population approach. We found that imipenem concentrations were best described by a two-compartment model in accordance with previously published studies [START_REF] Belzberg | Imipenem levels are not predictable in the critically ill patient[END_REF][START_REF] Novelli | Pharmacokinetic evaluation of meropenem and imipenem in critically ill patients with sepsis[END_REF][START_REF] Sakka | Population pharmacokinetics and pharmacodynamics of continuous versus short-term infusion of imipenem-cilastatin in critically ill patients in a randomized, controlled trial[END_REF]. The strength of our study is the number of patients with prospective collection of kinetic profiles with 6 points and a central laboratory for concentration assessment. To the best of our knowledge, the present study is the first study of imipenem PK using a population approach in the ICU.

In the population PK study performed by Lee et al. [START_REF] Lee | Comparison of 30-min and 3-h infusion regimens for imipenem/cilastatin and for meropenem evaluated by Monte Carlo simulation[END_REF] in healthy subjects, estimated imipenem clearance was 12.1 L/h and estimated central volume was 9.7 L. In ICU patients, we found a very similar clearance. The volume of distribution was estimated to 20.4 L in the final model, which is twice higher than that described in a healthy population. This increase is consistent with the clinical status of ICU patients. Indeed, inflammatory response in sepsis lead to an increased capillary For covariate selection we used a standard stepwise approach rather than a more modern approaches (such as the lasso method associated with cross-validation [START_REF] Ribbing | The lasso--a novel method for predictive covariate model building in nonlinear mixed effects models[END_REF]). Results of the selection steps were very consistent with our rich pharmacokinetic design. Of the 12 covariates studied, we found that CrCl4h, total body weight and serum albumin have a significant impact on the PK variability of imipenem. The addition of these three covariates reduces the variability of imipenem clearance and central Vd with a decrease of 10% for CL and 17% for V1. A recent study by Yoshizawa et al. (38) also showed the impact of creatinine clearance on imipenem clearance in patients with altered renal function. The other covariates tested, namely age and body weight, were not kept in their model.

Estimates of pharmacokinetic parameters