Ryan P C De Souza

Marcos V Moreira
email: moreira.mv@poli.ufrj.br

Jean-Jacques Lesage
email: jean-jacques.lesage@ens-paris-saclay.fr

A timed model for discrete event system identification and fault detection

Keywords: Discrete-event systems, System identification, Fault detection, Finite-state automata

We present in this paper a timed discrete event model for system identification with the aim of fault detection, called Timed Automaton with Outputs and Conditional Transitions (TAOCT). The TAOCT is an extension of a recent untimed model proposed in the literature, called Deterministic Automaton with Outputs and Conditional Transitions (DAOCT). Differently from the DAOCT, where only the logical behavior of the discrete event system is considered, the TAOCT takes into account information about the time that the events are observed, and, for this reason, it can be used for the detection of faults that cannot be detected by using untimed models, such as faults that lead the fault detector to deadlocks. The TAOCT represents the fault-free system behavior, and a fault is detected when the observed behavior is different from the behavior predicted by the model, considering both logical and timing information. A practical example is used to illustrate the results of the paper.

INTRODUCTION

In the past years, interest in the domain of fault diagnosis of discrete-event systems (DES) has increased. Since the introduction of the concept of diagnosability of discrete-event systems in [START_REF] Sampath | Diagnosability of discrete event systems[END_REF], several methods have been proposed in the literature for fault diagnosis of DES [START_REF] Debouk | Coordinated decentralized protocols for failure diagnosis of discrete event systems[END_REF][START_REF] Moreira | Polynomial time verification of decentralized diagnosability of discrete event systems[END_REF][START_REF] Zaytoon | Overview of fault diagnosis methods for discrete event systems[END_REF]Cabral et al., 2015). These methods provide a theoretical framework for the study of fault diagnosis of DES. However, their application to complex real systems, such as industrial plants, is a difficult task since these methods rely on an accurate model of the system, including its post-fault behavior.

In general, the modeling of a real-world system is very laborious and time consuming, since, depending on its size and complexity, it is very difficult to take into account all possible behaviors of the system in the model. This problem is amplified when we consider the post-fault behavior, since there may exist unpredicted consequences to a fault occurrence. In addition, the modeling process requires engineers that know the plant behavior, and are familiar with discrete-event modeling techniques. All these problems restrict the application of methods based on the complete system model to small systems, where the faultfree and faulty behaviors can be completely known.

In order to circumvent the aforementioned problems, fault detection techniques based on an identified fault-free model of the system have been proposed in the literature [START_REF] Roth | The concept of residuals for fault localization in discrete event systems[END_REF][START_REF] Klein | Fault detection of discrete event systems using an identification approach[END_REF]Moreira and Lesage, 2019a). In these works, the two main ideas are: (i) to automate the process of obtaining the fault-free model of the system by using identification;

This study was financed in part by CNPq, FAPERJ, and the Coordenac ¸ão de Aperfeic ¸oamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001. and (ii) when a fault has been detected through a discrepancy between the system behavior and the model, to use a technique based on residuals for fault localization.

A model for the identification of closed-loop industrial DES, called Non-Deterministic Autonomous Automaton with Outputs (NDAAO), is presented in [START_REF] Klein | Fault detection of discrete event systems using an identification approach[END_REF]. The identification procedure is based on the acquisition of binary signals exchanged between the programmable logic controller (PLC) and the plant. These signals correspond to sensor readings (inputs to the controller) and actuator commands (outputs of the controller), as shown in Fig. 1. Using the recorded data, an automaton which is capable of simulating the observed behavior is constructed. In [START_REF] Klein | Fault detection of discrete event systems using an identification approach[END_REF], with a view to obtaining a compact model, loops are introduced in the identified model, leading to the generation of sequences that have not been observed during the identification process. These sequences form the exceeding language generated by the identified model, and can be associated to non-detectable faults, since they are included in the fault-free model.

In order to reduce the exceeding language while keeping the compactness of the identified model, in Moreira and Lesage (2019a) a new automaton model, called Deterministic Automaton with Outputs and Conditional Transitions (DAOCT), is proposed. Conditions for the transposition of the transitions associated with the observed paths are added to the model by using a path estimation function. The use of the path estimation function reduces the exceeding language generated by the model in comparison with the NDAAO, and, consequently, it reduces the number of non-detectable faults. In [START_REF] Moreira | Fault diagnosis based on identified discrete-event models[END_REF], an algorithm for fault diagnosis using the DAOCT model proposed in Moreira and Lesage (2019a) It is important to remark that some faults cannot be detected by using untimed models. For instance, faults that prevent the system from generating events that are expected by the fault detector, lead it to a deadlock that makes the fault non-detectable. Thus, in some cases, it is necessary to add timing information in order to improve the capability of detecting faults. In [START_REF] Schneider | Determination of timed transitions in identified discrete-event models for fault detection[END_REF], the NDAAO model is extended to consider timing information, and timing constraints in the form of guards are added to each transition. Each transition has a guard that is defined as a single timing interval, and the guard pre-condition is satisfied when the event occurs in this interval. This model, called Timed Autonomous Automaton with Outputs (TAAO), is capable of detecting faults if the system deadlocks.

In this paper we introduce a new timed model for DES identification with the aim of fault detection. The new model, called Timed Automaton with Outputs and Conditional Transitions (TAOCT), is an extension of the DAOCT model in which timing information is added to the transitions in the form of guards for the transposition of the transitions. Since the TAOCT is based on the DAOCT, it inherits the same advantages of the underlying DAOCT model in terms of exceeding language reduction in comparison with the TAAO model, which is based on the NDAAO. In addition, differently from the TAAO model, the procedure to define the guards of the transitions of the TAOCT is capable of distinguishing different timing intervals associated with the same transition. Thus, more accurate guards can be defined, increasing the number of detectable faults in comparison with the TAAO model that considers a single timing interval for each transition. In addition, the timing information can also be used to refine the path estimate carried out by the path estimation function of the underlying DAOCT, reducing the number of non-detectable faults. It is important to remark that, as in the TAAO, deadlocks caused by faults can also be detected by using the TAOCT.

This paper is structured as follows. In Section 2, we present some preliminary concepts. Then, in Section 3, we explain the identification procedure. The TAOCT is defined in Section 4.

A practical example is presented and discussed in Section 5. Finally, conclusions are drawn in Section 6.

PRELIMINARIES

Let G be a Timed Automaton with Timing Intervals (TATI) defined as follows [START_REF] Cassandras | Introduction to Discrete Event Systems[END_REF]. Definition 1. A timed automaton with timing intervals is a sixtuple: G = (X, Σ, f , c g , guard, x 0 ), where X is the set of states, Σ is the finite set of events, f : X × Σ * → X is the transition function, c g is the global clock with value c g (t) ∈ R + , t ∈ R + , guard : X × Σ → A is the guard function, where A is the set of admissible timing intervals for the global clock c g , and x 0 ∈ X is the initial state.

2

A timed automaton with timing intervals is an automaton to which a global clock c g is added. Thus, timed sequences are generated by the TATI. Function guard : X × Σ → A specifies the timing conditions that need to be satisfied on the global clock for the transition to occur, and A are the admissible clock constraints defined as timing intervals. In the TATI, the global clock is reset to zero each time an event occurs. Let t denote the time that a state x ∈ X is reached in the timed automaton, and let τ denote the time that has elapsed after reaching state x. Then, transition (x, σ , x ), where x = f (x, σ ), occurs in the TATI if event σ occurs at time t + τ, and τ ∈ guard(x, σ ).

The TATI is a simplified version of the timed automaton with guards [START_REF] Cassandras | Introduction to Discrete Event Systems[END_REF]. In this paper, a similar formalism is used to add timing information to the proposed identification model.

A timed path of a timed automaton with timing intervals G is a sequence of states, events and time values that can be executed by G. A timed path p can be written in the following form:

p = (x 1 , σ 1 , τ 1 , x 2 , σ 2 , τ 2 , ..., x l-1 , σ l-1 , τ l-1 , x l ) , (1
) where x j ∈ X, j = 1, ..., l, are the states visited in the path, σ j ∈ Σ are the events that trigger the transition from state x j to x j+1 , j = 1, ..., l -1, and τ j ∈ R + are the time durations where the system remains in state x j up to the occurrence of event σ j , j = 1, ..., l -1. The length of a path p, denoted by |p|, is equal to the number of vertices in the path, i.e., |p| = l.

The timed sequence associated with timed path

p = (x 1 , σ 1 , τ 1 , x 2 , σ 2 , τ 2 , ..., x l-1 , σ l-1 , τ l-1 , x l ) is the sequence formed of events σ and time durations τ, s t = (σ 1 , τ 1 )(σ 2 , τ 2 ) . . . (σ l-1 , τ l-1 ).
It is important to remark that, in general, in the literature, timed sequences consider an absolute time that specifies when the event has occurred with respect to the beginning of the sequence execution, and not with respect to the sojourn time of each state of the path. However, the same information is contained in both types of timed sequences, since it is always possible to obtain a timed sequence with respect to the beginning of the process, from the timed sequence considering time durations s t , and vice-versa.

The set of non-negative integers is denoted by N, and the set formed only with 0 and 1 is denoted by N 1 = {0, 1}. The symbol '!' denotes 'is defined'.

IDENTIFICATION PROCEDURE

Let us consider a closed-loop industrial DES whose scheme is shown in Fig. 1. Suppose the controller has n i input signals, which correspond to sensor readings, and n o output signals, which correspond to actuator commands. Let us define the I/O vector, whose elements are the current status of each one of the controller signals: Suppose that at a given time instant t j , j ∈ N, the I/O vector becomes u j := u(t j ). Then, this vector remains unchanged until at least one of the controller signals changes its value at time t j+1 , leading to I/O vector u j+1 = u(t j+1 ). The events of the identified model are defined as the observed changes in the controller signals. Thus, when the I/O vector changes from u j to u j+1 , we consider that an event σ j has occurred. The sojourn time τ j is the time duration during which the I/O vector remained equal to u j , i.e., τ j = t j+1t j . Thus, when the observed I/O vector changes from u j to u j+1 , we say that transition (u j , σ j , τ j , u j+1 ) has been observed. A finite sequence of such transitions constitutes an observed timed path executed by the system p

u(t) := [i 1 (t) i 2 (t) ... i n i (t) o 1 (t) o 2 (t) ... o n o (t)] T , where i m (t) ∈ N 1 , m ∈ {1, ..., n i },
= (u 1 , σ 1 , τ 1 , u 2 , σ 2 , τ 2 , ..., u l-1 , σ l-1 , τ l-1 , u l ).
It is assumed in this paper that the system has a unique initial state, whose corresponding I/O vector is denoted as u 0 , and that all observed paths start at the initial state of the system, with the same I/O vector u 0 . The initial state may correspond, for example, to the beginning of a production cycle in the case of an industrial process. The objective of system identification is to compute a model that simulates the observed behavior described by the paths generated by the system. When untimed models are used, the behavior of the system is described only by the sequences of events that it generates. However, in the case of timed models, the identified model must simulate the observed sequences of events and must also be coherent with the time these events occur, i.e., it must be capable of simulating the timed sequence of events of the system. In both cases, it is possible that the identified model generates more sequences than those observed. In this case, there is an exceeding language that is associated with non-detectable faults if their sequences cannot be generated by the original fault-free system. In order to obtain an accurate model for identification, its exceeding language must be reduced.

It is important to remark that, in order to obtain a model by identification of a DES, capable of reproducing all possible behaviors of the system in fault-free operation, it is necessary to observe the paths executed by the system for an arbitrarily long time. However, in a finite time, only part of the paths generated by the system can be observed. This fact has a direct impact on the accuracy of the identified model and, therefore, on fault detection, since an observed behavior that the model is not capable of executing is interpreted as being faulty. If it is a fault-free behavior, then the fault detection system generates a false alarm. Thus, in order to reduce the number of false alarms, as in Moreira and Lesage (2019a), we assume in this paper that the system has been observed for a sufficiently long time such that all untimed paths of length up to a given number n 0 ∈ N have been observed. Moreover, we assume that each untimed path has been observed a number of times sufficient to capture the timing information regarding event occurrences.

The set formed of all observed events during the identification procedure is denoted by Σ. The set of all I/O vectors that have been observed is denoted by Ω. Thus, Ω ⊆ N n i +n o 1 .

Time-interval paths

Let us assume that N timed paths (not necessarily distinct) have been observed. Each path is denoted by

p n = u n,1 , σ n,1 , τ n,1 , u n,2 , σ n,2 , τ n,2 , ..., u n,l n -1 , σ n,l n -1 , τ n,l n -1 , u n,l n ,
where n ∈ {1, ..., N}, u n, j ∈ Ω, for j = 1, ..., l n , and σ n, j ∈ Σ and τ j,n ∈ R + , for j = 1, ..., l n -1.

The set containing all of the N observed paths is denoted by P. Let us partition P as follows: P = P 1 ∪P 2 ∪... ∪P r , r ≤ N, where all paths that form P i share the same logical sequence of states and events, given by a unique untimed path p u i = (u i,1 , σ i,1 , u i,2 , σ i,2 , ..., u i,l i -1 , σ i,l i -1 , u i,l i ) for each set P i , i = 1, . . . , r. As in Moreira and Lesage (2019a), we assume here that none of the untimed paths p u i has an associated sequence of events s i = σ i,1 σ i,2 . . . σ i,l-1 that is a prefix of the sequence of events of another path p u j , where i = j.

For each set of paths P i we build a time-interval path p i given by:

p i = u i,1 , σ i,1 , I i,1 , ..., u i,l i -1 , σ i,l i -1 , I i,l i -1 , u i,l i ,
where u i, j = u n, j , j = 1, . . . , l i , σ i, j = σ n, j , j = 1, . . . , l i -1, for any n ∈ {1, . . . , N} such that p n ∈ P i , and I i, j ⊂ R + , for j = 1, . . . , l i -1. Note that, in the time-interval path p i , instead of timing values τ i, j , sets I i, j ⊂ R + are used. Each set I i, j is built in the following manner: (i) construct set T i, j = {τ n, j : p n ∈ P i };

(ii) form a set of clusters S = {S 1 , . . . , S K } from T i, j ; (iii) form set

I i, j = K h=1 [max{0, min S h -δ }, max S h + δ ].
Note that an uncertainty δ ∈ R + is added to form each real interval in I i, j to represent possible errors in measuring the time occurrence of the events, caused, for example, by the scan cycle of the controller. The value of δ must be chosen such that the intervals do not overlap with each other, i.e., the user-specified threshold for defining the separation of the clusters must be greater than 2 × δ .

The cluster function, used in the process of forming the sets I i, j , takes as input a finite set of real values, and gives as output a set of subsets S, each one containing the values which are more densely grouped compared with the values in the other subsets of the input set. We say that the values in each subset S h ∈ S form a cluster. Clustering methods have been widely discussed in the literature, and there are several methods that can be used to implement the cluster function in the process of obtaining the sets I i, j , as it can be seen in [START_REF] Jain | Data clustering: A review[END_REF]. In this paper, the NEAREST NEIGHBORS method [START_REF] Lu | A sentence-to-sentence clustering procedure for pattern analysis[END_REF] is used. In this method, an element is assigned to a cluster if the minimum distance between this element and the other members of the cluster is smaller than a user-specified threshold.

Parametric identification approach

In Moreira and Lesage (2019a), inspired by [START_REF] Klein | Fault detection of discrete event systems using an identification approach[END_REF], the parametric identification algorithm allows to obtain a model satisfying an important property called k-completeness that guarantees that, once the value of a free parameter k has been chosen for system identification, all and only all sequences of length 1 to k + 1 that have been observed are generated by the identified DAOCT. By increasing the value of the free parameter k, the exceeding language generated by the DAOCT reduces, but the size of the model grows. Thus, there is a trade-off to be found between complexity and accuracy of the identified model.

Since the timed model proposed in this paper is based on the DAOCT model presented in Moreira and Lesage (2019a), we use the same strategy proposed in Moreira and Lesage (2019a) to obtain a parameterized model for identification. In order to do so, the same steps for the computation of the DAOCT model are considered here, where the first step is the computation of modified paths p k i , according to the free parameter k, as follows:

p k i = y i,1 , σ i,1 , I i,1 , ..., y i,l i -1 , σ i,l i -1 , I i,l i -1 , y i,l i , where
y i, j = (u i, j-k+1 , ..., u i, j ), if k ≤ j ≤ l i (u i,1 , ..., u i, j ), if j < k .
Set Ω k is defined as Ω k := i∈R {y i, j : j = 1, . . . , l i }. The last element of y i, j is denoted by y l i, j . Thus, y l i, j = u i, j , j = 1, . . . , l i , and i ∈ R.

TIMED AUTOMATON WITH OUTPUTS AND CONDITIONAL TRANSITIONS

We present in the sequel the timed model proposed in this paper for the identification of DES. The model is based on Definition 1 of Timed Automaton with Timing Intervals, and on the DAOCT model proposed in Moreira and Lesage (2019a). Definition 2. The Timed Automaton with Outputs and Conditional Transitions is a ten-tuple: TAOCT = (X, Σ, f , c g , Ω, λ , R, g, x 0 , X f ), where X is the finite set of states, Σ is the finite set of events, f : X × Σ * → X is the transition function, c g is the global clock, with value c g (t

) ∈ R + , t ∈ R + , Ω ⊆ N n i +n o 1
is the set of outputs, where n i and n o denote, respectively, the number of inputs and outputs of the controller, λ : X → Ω is the state output function, R = {1, 2, ..., r} is the set of path indexes, g : X × Σ × R → C is the guard function, x 0 ∈ X is the initial state, and X f ⊆ X is the set of final states.

2

The set of admissible constraints C is formed of all sets I ⊂ R + . As in the TATI, presented in Definition 1, in the TAOCT a unique global clock c g is used, and it is reset every time a transition occurs. Function g(x, σ , i) specifies a subset of R + to which the clock value c g (t) must belong so that transition (x, σ , x ), where x = f (x, σ ), associated with observed path p i , can occur. The output function λ is the same presented in the DAOCT model, and associates an I/O vector with each state x ∈ X of the model.

Differently from the DAOCT, where a path estimation function is defined in the model, in the TAOCT the path estimation function θ : X × Σ × R + → 2 R can be defined using the guard function g as follows:

θ (x, σ , τ) = {i ∈ R : τ ∈ g(x, σ , i)}. Function θ (x, σ , τ) provides the path estimate when the current state is x, and event σ occurs at clock value c g (t) = τ. If f (x, σ ) is not defined for a given path j, then g(x, σ , j) is not defined, and j ∈ θ (x, σ , τ) for any value of τ. Moreover, if f (x, σ ) is defined, but the observed time τ does not belong to g(x, σ , i), then i ∈ θ (x, σ , τ).

Algorithm 1 presents the procedure for constructing the identified TAOCT model from the modified time-interval paths p k i , for i = 1, . . . , r. In order to do so, function λ : X → Ω k is defined. The guards and its associated path estimation function are used in the TAOCT evolution rule to reduce the exceeding language generated by the model, reducing the number of non-detectable faults. We present in the sequel the timed language generated by the model. In order to do so, we first define the set of timed events Σ t of the TAOCT as follows:

Σ t := (σ , τ) ∈ Σ × R + : (∃x ∈ X)[θ (x, σ , τ) = / 0] . Set Σ t is
formed of all pairs (σ , τ), where σ ∈ Σ and τ ∈ R + , such that there exist x ∈ X, where f (x, σ ) is defined, and at least one path p i , such that τ ∈ g(x, σ , i). Let us define Σ * t as the set of all possible sequences formed of elements σ t ∈ Σ t , and assume that the empty sequence ε belongs to Σ * t . Function ψ : Σ * t → Σ * removes the timing information from a timed sequence in Σ * t , obtaining its equivalent untimed se-Algorithm 1: TAOCT identification

Input: Modified time-interval paths p k i , for i = 1, ..., r Output: TAOCT = (X, Σ, f , c g , Ω, λ , R, g, x 0 , X f ) 1 Create initial state x 0 and define λ (x 0 ) = λ (x 0 ) = y 1,1 2 X ← {x 0 }, Σ ← / 0, Ω ← {y 1,1 }, R ← {1, ..., r}, X f ← / 0 3 Create global clock c g 4 for i = 1 to r do 5 for j = 1 to l i -1 do 6 Find state x ∈ X such that λ (x) = y i, j 7 if λ (x) = y i, j+1 ∀x ∈ X then 8 Create state x and define λ (x ) = y i, j+1 9 X ← X ∪ {x } 10 Σ ← Σ ∪ σ i, j 11 λ (x ) ← y l i, j+1 12 Ω ← Ω ∪ {λ (x )} 13 else 14 Find x ∈ X such that λ (x ) = y i, j+1 15 f (x, σ i, j ) ← x 16 if g(x, σ i, j , i)! then 17 g(x, σ i, j , i) ← g(x, σ i, j , i) ∪ I i, j 18 else 19 g(x, σ i, j , i) ← I i, j 20 if j = l i -1 then 21 X f ← X f ∪ {x } quence, i.e., for any s t = (σ 1 , τ 1 )...(σ n , τ n ) ∈ Σ * t , then ψ(s t ) = σ 1 ...σ n ∈ Σ * . By convention, ψ(ε) := ε.
Let us define, now, recursively, the extended path estimation function θ s : X × Σ * t → 2 R as follows: θ s (x, ε) = R, and for any sequence s t (σ , τ) ∈ Σ * t , where s t ∈ Σ * t and (σ , τ) ∈ Σ t , we have that:

θ s (x, s t (σ , τ)) =    θ s (x, s t ) ∩ θ (x , σ , τ), where x = f (x, ψ(s t )), if f (x, ψ(s t )σ )! undefined, otherwise.
The timed language generated by the TAOCT model is given by:

L t := {s t ∈ Σ * t : θ s (x 0 , s t ) = / 0} .
Language L t simulates the timed sequences of the paths used in the identification procedure as stated in the sequel. Theorem 1. Let L t,Obs denote the language formed of all timed sequences associated with the observed paths p n , for n = 1, 2, . . . , N. Then, L t,Obs ⊆ L t . 2

Proof. The proof is straightforward from the construction of the TAOCT according to Algorithm 1, and will be omitted due to the lack of space.

Since timing information is added to the TAOCT model, the fault detection scheme based on the timed model can detect faults that are not possible to be detected by using the DAOCT model. There are three faulty scenarios for which the fault can be detected by the TAOCT and not by the DAOCT model: (i) faults that lead the system to a deadlock; (ii) faults that cause the occurrence of an event σ in a state x, such that f (x, σ ) is defined, at a time instant τ that does not belong to any set defined by the guard conditions g(x, σ , i), for all i ∈ R; and (iii) faults that cause the occurrence of an event σ at a time instant τ that satisfies a guard condition, but leads the path estimation function to θ s (x 0 , s t ) = / 0. The first type of fault can be detected by counting the time elapsed after reaching a state of the model. If the elapsed time is greater than all possible times defined in the guards, then the deadlock is detected. The second and third types of faults can be detected by verifying if the observed timed sequence can be generated by the model, since, in both cases, it can be seen that, although f (x 0 , ψ(s t )) is defined, we have that θ s (x 0 , s t ) = / 0, i.e., the observed timed sequence is not in L t .

PRACTICAL EXAMPLE

The identification method proposed in this paper is illustrated using the sorting unit system presented in Figure 2. Three different types of pieces are sorted in the system: white plastic pieces (WP), black plastic pieces (BP), and metallic pieces (M). Each type of piece is pushed to one of the three slides shown on the bottom of Figure 2, such that pieces of type WP are pushed to the right slide, pieces of type M are pushed to the slide in the middle, and pieces of type BP are pushed to the left slide.

On the right of Figure 2, there is a stack magazine where the pieces are stored in any order. In the sorting process, the pieces at the bottom of the stack magazine are placed onto the conveyor belt by a pneumatic pusher. Then, the conveyor belt is turned on, and the piece is moved in the direction of two sensors in order to determine its type. An inductive sensor detects metallic pieces (type M), and an optical sensor detects metallic (type M) and white plastic pieces (type WP). If a black plastic piece (type BP) is on the conveyor, then none of the sensors is capable of detecting it. The optical sensor is located close to the inductive sensor, such that metallic pieces are detected by both sensors almost at the same time.

It is also important to remark that there is a photoelectric sensor next to each sorting pusher on the conveyor. When a piece is detected by the photoelectric sensor next to the pusher that should remove it from the conveyor, the conveyor is stopped and the pusher is extended. Then, the pusher is retracted and a new piece can be placed on the conveyor by the pusher of the stack magazine. The sorting unit system has 13 sensors and 6 actuator signals.

Thus, the controller has 19 input and output signals. The initial state of all observed paths is defined as the I/O vector corresponding to the case where the conveyor belt is turned off, and all pushers are retracted.

During the identification process, 2294 timed paths were observed, which corresponds to two hours and fifty three minutes of observation of the controller signals. From all the observed timed paths, 12 logically different paths were obtained. As there are only three types of pieces, we could expect only three different logically distinguishable sequences. However, it has been observed that, since the inductive and optical sensors are very close to each other, the order of sensor readings (rising and falling edges) may change for different sorting cycles of metallic pieces, increasing the number of paths.

In this practical example, we have chosen the threshold for the clustering algorithm to be 100ms, since this value is coherent with the dynamics of the sorting system. In this case, the cluster function returned a single closed interval for each transition of the time interval paths p i , for i = 1, . . . , 12.

After obtaining the time interval paths p i , the TAOCT model is identified following the steps of Algorithm 1. In this case, for k = 1, the identified TAOCT has 26 states and 38 transitions.

In Figure 3, we present only part of the identified TAOCT due to lack of space. Since each transition of p i is associated with a single time interval, then all guards g(x, σ , i) of the identified TAOCT are single closed intervals. In Figure 3 we do not present the guards due to lack of space. The events of the TAOCT are rising and falling edges of the elements of the I/O controller vector, which are represented by ↑z and ↓z, respectively, where z is the position of the signal in the I/O vector. It is important to remark that an event can be formed of more than one rising or falling edge of controller signals.

In the sequel, we illustrate the three faulty scenarios for which the fault can be detected thanks to the timing information added to the TAOCT model.

In the first scenario, let us consider that, after a piece is placed on the conveyor belt by the pusher of the stack magazine, the pusher stuck extended and cannot be retracted. The path associated with this behavior is p = (x 0 , ↑17, x 1 , ↓2, x 2 , ↑1. ↓17. ↑18, x 3 ), where 17, 2, 1, and 18 are, respectively, the command to extend the pusher of the stack magazine, the sensor that indicates that the pusher of the stack magazine is retracted, the sensor that indicates that the pusher of the stack magazine is extended, and the command to retract the pusher. In this case, the system deadlocks, since the conveyor belt is turned on only after the retraction of the pusher is detected (↑ 2). This fault cannot be detected by using the DAOCT model, but can be detected by using the TAOCT model, since the falling edge of the sensor that indicates that the pusher is extended (↓1) must occur before the maximum time of the guard g(x 3 , ↓1, 1) = [87, 161]. Thus, when the elapsed time is greater than 161 milliseconds, the fault is detected.

To illustrate the second faulty scenario, let us consider a fault in the speed controller of the conveyor that makes it work faster than expected. Let us consider, for instance, that after a BP piece is placed on the conveyor belt, and the conveyor is turned on (↑19), which corresponds to state x 5 of the TAOCT, it reaches the photoelectric sensor next to the first sorting pusher (↑13) in a time smaller than the minimum time of the guard g(x 5 , ↑13, 4) = [4058,4146]. Thus, the fault is detected. It is important to remark that, since ↑13 is coherent with the logical behavior of the system, then the diagnosis system based on the DAOCT model would not be capable of detecting the fault.

The third faulty scenario can be illustrated by the following example. Consider all observed paths associated with BP or M pieces, and consider that the piece is in front of the photoelectric sensor next to the right sorting pusher, i.e., the rising edge of the photoelectric sensor (↑ 13) has been observed, which corresponds to state x 8 of Figure 3. By analyzing the time elapsed between ↑13 and ↓13, two distinct sets of time values can be defined according to the type of piece, as shown in Figure 4. This occurs because metallic pieces are detected for a longer time than plastic pieces by the photoelectric sensor due to their brightness. Let us consider now a fault that causes both the optical and inductive sensors to fail at the same time. In this case, a piece of type M would lead to the same logical behavior than a BP piece, making such a fault nondetectable by the diagnosis system based on the DAOCT model. However, as shown in Figure 4, it is possible to distinguish the types of pieces by using the guards associated with the timed paths. While a BP piece would take some time in the interval g(x 8 , ↓13, 4) = [915, 937], a metallic piece would take some value in the interval [1035, 1052] milliseconds, which corresponds to the union of all guards defined in state x 8 , for event ↓13, and the timed paths associated with metallic pieces. Thus, if a metallic piece is on the conveyor belt, and the fault occurs, the time elapsed between the rising and falling edges of the photoelectric sensor will be compatible with the guard condition associated with metallic pieces, and non-compatible with the guard condition of black plastic pieces. Since the logical behavior is not coherent with the timing behavior, the fault is detected.

CONCLUSIONS

In this paper, a new timed model for the identification of DES with the aim of fault detection is proposed. In this model, called Timed Automaton with Outputs and Conditional Transitions (TAOCT), timing information regarding the occurrence of events are added as guards to the transitions. By doing so, a refinement of the path estimation can be carried out using the timing information, and faults that cannot be detected by using untimed models, can be detected by using the TAOCT. We are currently working on a fault detection and isolation algorithm based on the TAOCT.

Cabral, F.G., Moreira, M.V., Diene, O., and Basilio, J.C. (2015). A Petri net diagnoser for discrete event systems modeled

  Fig. 1. Closed-loop system showing the signals exchanged between plant and controller.

  are the values of the input signals at time t ∈ R + , and o m (t) ∈ N 1 , m ∈ {1, ..., n o }, are the values of the output signals at time t ∈ R + .

Fig. 2 .

 2 Fig. 2. Sorting unit system of the practical example.

Fig. 3 .

 3 Fig. 3. Part of the TAOCT model obtained for the practical example.

Fig. 4 .

 4 Fig. 4. Histogram showing the observed times elapsed between the rising edge and falling edge of the photoelectric sensor next to the first sorting pusher (BP in black and M in grey) in the observed timed paths.


Mélanie Prague

Linda Wittkop

Quentin Clairon

Dan Dutartre

Rodolphe Thiébaut
email: rodolphe.thiebaut@u-bordeaux.fr

Boris P Hejblum

Population modeling of early COVID-19 epidemic dynamics in French regions and estimation of the lockdown impact on infection rate

Keywords: Compartmental model, COVID-19, France, Population modeling, Non-pharmaceutical intervention, SARS-CoV-2 Declarations of interest: none Highlights

come   L'archive ouverte pluridisciplinaire

Introduction

In December 2019, grouped pneumonia cases have been described in the Hubei province, China and SARS-CoV2 was identified on January, 7 th as the cause of this outbreak (Li et al., 2020a;[START_REF] Zhu | A Novel Coronavirus from Patients with Pneumonia in China, 2019[END_REF]. SARS-CoV2 causes the viral disease which has been named COVID-19 (World Health Organization, 2020b).

SARS-CoV2 rapidly spread all over the world and the pandemic stage was declared on March 11 th by the World Health Organization (2020c). On April 28 th , over 1,773,084 cases (in accordance with the applied case definitions and testing strategies in the affected countries) including 111,652 deaths have been reported (World Health Organization, 2020a).

The first case in France was declared on January, 24 th (Bernard-Stoecklin et al., 2020) and on April 13 th , Santé Publique France reported 98,076 confirmed cases and 14,967 hospital deaths due to COVID-19.

COVID-19 includes non-specific symptoms such as fever, cough, headache, and specific symptoms such as loss of smell and taste [START_REF] Gane | Isolated sudden onset anosmia in covid-19 infection. a novel syndrome?[END_REF][START_REF] Greenhalgh | Covid-19: a remote assessment in primary care[END_REF]. The virus is transmitted through droplets and close unprotected contact with infected cases. The majority (around 80 %) of infected cases have a mild form (upper respiratory infection symptoms) without specific needs in terms of care. Around 20 % of cases need hospitalization and among those are severe forms (severe respiratory distress) which will need to be admitted to intensive care units (ICU) with potential need of mechanical ventilation. The percentage of patients in need for ICU care varies between 5 % reported from China [START_REF] Wj | Clinical characteristics of coronavirus disease 2019 in china[END_REF] and 16 % reported from Italy [START_REF] Grasselli | Critical Care Utilization for the COVID-19 Outbreak in Lombardy, Italy: Early Experience and Forecast During an Emergency Response[END_REF]. The number of ICU beds in France was 5,058 at the end of 2018 (DREES, 2019) (although it is currently being increased, having doubled and aiming to reach 14,000 according to the French minister of Health). Thus, the availability of ICU beds with mechanical ventilation is one of the major issues as facilities are not prepared to deal with the potential increase of the number of patients due to this epidemic.

Unprecedented public-health interventions have been taken all over the world [START_REF] Kraemer | The effect of human mobility and control measures on the covid-19 epidemic in china[END_REF] to tackle this epidemic. In France, interventions such as heightening surveillance with rapid identification of cases, isolation, contact tracing, and follow-up of potential contacts were initially implemented. But as the epidemic continued growing, comprehensive physical distancing measures have been applied since March 15 th , 2020 including closing of restaurants, non-vital business, schools and universities etc, quickly followed by state-wide lockdown on March 17 th 2020. The president has announced on April 13 th 2020, a progressive lifting of the lockdown from May 11 th 2020 onwards. In Wuhan (Hubei, China), the extremely comprehensive physical distancing measures in place since January 23 rd have started to be relaxed after 2 months of quarantine and lifted completely on April 8 th 2020 [START_REF] Tian | An investigation of transmission control measures during the first 50 days of the COVID-19 epidemic in china[END_REF][START_REF] Wu | Characteristics of and Important Lessons From the Coronavirus Disease 2019 (COVID-19) Outbreak in China: Summary of a Report of 72,314 Cases From the Chinese Center for Disease Control and Prevention[END_REF].

Interestingly, these interventions have been informed by mathematical models used to estimate the epidemic key parameters as well as unmeasured compartments such as the number of infected people. Another interesting outcome is the forecast of the COVID-19 epidemic according to potential interventions.

Several models have already been proposed to model and forecast the COVID-19 epidemic using compartment models [START_REF] Fang | Transmission dynamics of the covid-19 outbreak and effectiveness of government interventions: A data-driven analysis[END_REF][START_REF] Tang | The effectiveness of quarantine and isolation determine the trend of the COVID-19 epidemics in the final phase of the current outbreak in china[END_REF][START_REF] Wang | Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan[END_REF] or agent based models [START_REF] Domenico | Report #8: Expected impact of school closure and telework to mitigate covid-19 epidemic in france[END_REF][START_REF] Ferguson | Report 9: Impact of non-pharmaceutical interventions (npis) to reduce covid19 mortality and healthcare demand[END_REF][START_REF] Wilder | The role of age distribution and family structure on covid-19 dynamics: A preliminary modeling assessment for hubei and lombardy (march 31[END_REF], its potential impact on intensive care systems [START_REF] Fox | Modelling the impact of COVID-19 upon intensive care services in New South Wales[END_REF][START_REF] Massonnaud | Covid-19: Forecasting short term hospital needs in france[END_REF], and to estimate the effect of containment measurements on the dynamics of the epidemic [START_REF] Magal | Predicting the number of reported and unreported cases for the covid-19 epidemic in south korea, italy, france and germany[END_REF][START_REF] Prem | The effect of control strategies to reduce social mixing on outcomes of the covid-19 epidemic in wuhan, china: a modelling study[END_REF]. Most of those rely on simulations with fixed parameters and do not perform direct statistical estimations from incident data [START_REF] Massonnaud | Covid-19: Forecasting short term hospital needs in france[END_REF]. [START_REF] Roques | Using early data to estimate the actual infection fatality ratio from COVID-19 in france[END_REF] used French national data but did not use a population approach to model the epidemic at a finer geographical granularity. Yet, the dynamics of the epidemic can be very heterogeneous between regions inside a given country resulting in tremendous differences in terms of needs for hospital and ICU beds [START_REF] Massonnaud | Covid-19: Forecasting short term hospital needs in france[END_REF]. Moreover, the data collection yields noisy observations, that we deal with a statistical modeling of the observation process, rather than altering the data by e.g. smoothing such as in [START_REF] Roques | Using early data to estimate the actual infection fatality ratio from COVID-19 in france[END_REF].

In the present study, we use public data from the COVID-19 outbreak in France to estimate the dynamics of the COVID-19 epidemic in France at the regional level. We model the epidemic with a SEIRAH model, which is an extended Susceptible-Exposed-Infectious-Recovered (SEIR) model accounting for time-varying population movements, non-reported infectious subjects (A for unascertained) and hospitalized subjects (H) as proposed by [START_REF] Wang | Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan[END_REF] to model the epidemic in Wuhan. Parameters from the model are estimated at the regional scale using a population approach which allows for borrowing information across regions, increasing the amount of data and thereby strengthening the inference while allowing for local disparities in the epidemic dynamics. Furthermore, we use forward simulations to predict the effect of non-pharmaceutical interventions (NPI) (such as lift of lockdown) on ICU bed availability and on the evolution of the epidemic. Section 2 introduces the data, the model and the necessary statistical tools, Section 3 presents our results and Section 4 discusses our findings and their limits.

Methods

Because epidemics spread through direct contacts, their dynamics have a strong spatial component. While traditional compartment models do not account for spatiality, we propose to take it into account by: i) modeling the epidemic at a finer, more homogeneous geographical scale (this is particularly important once lockdown is in place); ii) by using a population approach with random effects across French regions which allows each region to have relatively different dynamics while taking all information into account for the estimation of model parameters iii) aligning the initial starting time of the epidemic for all regions. The starting date in each region was defined as the first date with incident confirmed cases of COVID-19 directly followed by 3 additional consecutive days with incident confirmed cases as well. This criterion of 4 consecutive days with incident cases is needed in particular for the Île-de-France region which had 3 consecutive days with 1 imported confirmed case in late January which did not lead to a spreading outbreak at that time.

Data sources

Open-data regarding the French COVID-19 epidemic is currently scarce, as the epidemic is still unfolding. Santé Publique France (SPF) in coordination with the French regional health agencies (Agences Régionales de Santé -ARS) has been reporting a number of aggregated statistics at various geographical resolutions since the beginning of the epidemic. During the first weeks of the epidemic in France, SPF was reporting the cumulative number of confirmed COVID-19 cases with a positive PCR test. Other French surveillance resources such as the Réseau Sentinelles [START_REF] Valleron | A computer network for the surveillance of communicable diseases: the french experiment[END_REF] or the SurSaUD R database [START_REF] Caserio-Schönemann | The French syndromic surveillance system SurSaUD R[END_REF] quickly shifted their focus towards COVID-19, leveraging existing tools to monitor the ongoing epidemic in real time, making as much data available as possible (given privacy concerns). In this study, we combined data from three different opendata sources: i) the daily release from SPF; ii) the SurSaUD R database that started recording visits to the Emergency room for suspicion of COVID-19 on February, 24 th ; iii) the Réseau Sentinelles which started estimating the weekly incidence of COVID-19 in each French region on March 16 th . From the daily release of SPF, we computed the daily incident number of confirmed COVID-19 cases (i.e. with a positive PCR test) in each region. In addition we used the incident number of visits to the emergency room for suspicion of COVID-19 in each region from the SurSaUD R database using the OSCOUR R network that encompasses more than 86% of all French emergency services [START_REF] Caserio-Schönemann | The French syndromic surveillance system SurSaUD R[END_REF]. Although this does not represent the full extend of hospitalized COVID-19 cases, it is the only public data available that early in the epidemic, when the majority of COVID-19 cases at hospitals were admitted through emergency rooms. Finally, we used the Réseau Sentinelles network's weekly incidence estimates of symptomatic cases (including non confirmed cases) to set the ratio between ascertained and unascertained cases in each region (later denoted as r i ). Table 1 presents these observed data. Of note, we studied the epidemic in the 12 Metropolitan French regions -excluding the Corsican region (Corse) which exhibits different epidemic dynamics, possibly due to its insular nature. [START_REF] Wang | Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan[END_REF] extended the classic SEIR model to differentiate between different statuses for infected individuals: ascertained cases, unascertained cases and cases quarantined by hospitalization. The model, assuming no population movement, is presented in Figure 1. The population is divided into 6 compartments: susceptible S, latent E, ascertained infectious I, unascertained infectious A, hospitalized infectious H, and removed R (recovered and deceased). This model assumes that infections are well-mixed throughout the population, ignoring any spatial structure or compartmentalization by population descriptors such as age. Such assumptions make it particularly relevant to infer the dynamics of the French epidemic at the regional level (a finer geographical scale at which such hypotheses are more likely to hold). Figure 1 illustrates the dynamics between those 6 compartments that are characterized by the following system of six Ordinary Differential Equations (ODE):

Model

Structural model of the epidemic

                                         dS dt = - bS(I + αA) N dE dt = bS(I + αA) N - E D e dI dt = rE D e - I D q - I D I dR dt = (I + A) D I + H D h dA dt = (1 -r)E D e - A D I dH dt = I D q - H D h (1)
Model parameters are described in Table 2. Of note, given a combination of parameters and initial states of the system ξ = (b, D q , r, α, D e , D I , D h , N, S(t = 0), E(t = 0), I(t = 0), R(t = 0), A(t = 0), H(t = 0)), using a solver of differential equations, it is possible to deterministically compute at any time t the quantities S(t, ξ), E(t, ξ), I(t, ξ), R(t, ξ), A(t, ξ), and H(t, ξ).

Observation model

Observation processes In our case, none of the compartments of the system are directly observed: the only observations considered are i) the number of daily incident infectious ascertained cases denoted Y 1 , and ii) the number of daily incident hospitalized infectious cases denoted Y 2 . These observations are the only one available both before and after the initiation of lockdown. Those two quantities are modeled in Equation ( 1) respectively as observations from the

I (in) (t, ξ) = rE(t,ξ)
De and H (in) (t, ξ) = I(t,ξ) Dq random variables, which are the numbers of new incident cases at time t given the parameters ξ in compartment I and H respectively. Because these are count processes, we propose to model their observations Y 1 and Y 2 with Poisson likelihoods:

P (Y 1 i = k 1 ) = e -I (in) (t,ξ) I (in) k 1 (t, ξ) k 1 ! (2) P (Y 2 i = k 2 ) = e -H (in) (t,ξ) H (in) k 2 (t, ξ) k 2 !
where k 1 and k 2 are the respective numbers of cases.

Initial values

The initial states of all compartments at the date of epidemic start (t = 0) for region i are also important drivers of the dynamics. Some of them can be approximated by quantities directly depending on the observations: i)

I i (t = 0) = t≤0 Y 1 i (t), ii) H i (t = 0) = t≤0 Y 2 i (t)
, and iii) R i (t = 0) = 0. Others, namely compartments A and E are not directly observed, and we evaluate these initial quantities. Due to variation in data collection protocols and the initial size of regional outbreaks, this estimation is particularly important. Indeed, the number of daily incident cases at t = 0 ranges from 1 to 37 cases depending on the region. A i (t = 0) is set as

I i (t=0)(1-r)
r by model assumption. E i (t = 0) is estimated as a direct parameter of the model (see Section 2.2.3). Finally, because the total population size is equal to N i , we have

S i (t = 0) = N i -E i (t = 0) -I i (t = 0) -R i (t = 0) -A i (t = 0) -H i (t = 0).

Statistical population model

The goal of this study is to model the epidemic of COVID-19 in France, but at the regional level using a population approach. This is done using a mixed effect model. In this inference framework, baseline parameters governing the dynamics of the epidemic in each region are assumed to be drawn from a shared distribution which allows for heterogeneity between regions, known as the random effects. We use the log-normal distribution for all parameters to ensure their positivity during estimation. Because public health policies changed over the time period of observation of the epidemic, we incorporate explanatory covariates such as physical distancing by lockdown (C 1 and C 2 ) as a time-dependent effect on the transmission of the disease b. Covariate C 1 is 0 until 2020-03-17 date of the start of the policy in France and then set to 1. Covariate C 2 is 0 until 2020-03-25 assuming that social distancing behaviours build up in a week. In other words, we have ∀i = 1, . . . , 12 (where i is the region identifier):

log(b i (t)) = b 0 + β 1 C 1 i (t) + β 2 C 2 i (t) + u b i , u b i ∼ N (0, σ 2 b ) log(D q i ) = D q 0 + u Dq i , u Dq i ∼ N (0, σ 2 Dq ) log(E i (t = 0)) = E(t = 0) + u E 0 i , u E 0 i ∼ N (0, σ 2 E 0 ) (3)
The parameters (b 0 , D q 0 , E(t = 0)) are mean shared values in the population, and can be seen as the country values for these parameters. The inter-region random-effect (u b i , u Dq i , u E 0 i ) are normally distributed and assumed independent. So the vector of parameters in the model for each region is

ξ i = (b i , D q i , β 1 , β 2 , r i , α, D e , D I , D h , N i , S i (t = 0), E i (t = 0), I i (t = 0), R i (t = 0), A i (t = 0), H i (t = 0)). β 1 can be interpreted through K 1 = exp(-β 1 )
which is the factor by which transmission is modified after the start of lockdown during the first week. The factor by which transmission is modified after that first week of confinement is given by K 2 = exp(-β 1 -β 2 ). The coefficients β = (β 1 , β 2 ) are expected to be negative as lockdown aims at reducing transmission. Interestingly, with our approach, we can evaluate whether or not there is a statistically significant effect of lockdown on the transmission by testing β = 0 using a Wald test.

Inference

Estimation

Region-specific model parameters Based on the results from the theoretical identifiability analysis of the structural model of the epidemic from Equation (1) (see Supplementary Materials S2), we estimate the parameters (b i , D q i , r i ) as well as the initial state (E 0i ) when the epidemic begins being reported in each region i. We used the Monolix software version 2019R2 (Lixoft SAS, 2019) to estimate those five parameters by maximizing the likelihood of the data given the model and the other fixed parameters. This software relies on a frequentist version of the Stochastic Approximation Expectation Maximization (SAEM) algorithm [START_REF] Delyon | Convergence of a Stochastic Approximation Version of the EM Algorithm[END_REF] and standard errors are calculated via estimation of the Fisher Information Matrix [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF], which is derived from the second derivative of the log-likelihood evaluated by importance sampling. In addition, we use profile-likelihood to confirm that no further information can be gained from the data at hand on parameters α, D e and D I by running the SAEM algorithm multiple times while setting these parameters to different values and obtaining similar maximum likelihood values (meaning more data would be needed to be able to estimate those parameters). During inference, practical identifiability of the model is evaluated by the ratio of the minimum and maximum eigenvalues of the Fisher Information Matrix, that will be referred as "convergence ratio" in the reminder of the manuscript. Convergence of the SAEM algorithm was assessed by running multiple SAEM chains and checking that they all mix around similar probability distributions.

Compartments dynamics

We are particularly interested in the trajectories of the model compartments. We use Monte Carlo methods (parametric bootstrapping) to compute the confidence intervals accounting for the uncertainty in estimating the structural and statistical model parameters. For all compartments C(t, ξ i ) (C being S, E, I, R, A, or H) the 95% confidence interval is estimated by sampling from the posterior distributions of the model parameters to simulate 1,000 trajectories, and taking the 2.5% and 97.5% percentiles of these simulated trajectories. We also added to it the error measurement given by the Poisson distribution of the outcomes. Other outcomes of interest are the number of ICU beds needed and the number of death (D) in a given region at a given time. These quantities are not specifically modeled by our mathematical structural model. However, it is possible to roughly approximate them by assuming that they represent a percentage of the hospitalized cases H(t, ξ i ) and removed cases R(t, ξ i ). We assume that ICU (t, ξ i ) = 0.25 × H(t, ξ i ) which is consistent with the prevalence of ICU cases among hospitalized cases at the French national level. Based on the estimation of the Infection Fatality Ratio (IFR) from [START_REF] Roques | Using early data to estimate the actual infection fatality ratio from COVID-19 in france[END_REF], we get a rough estimation of D(t, ξ i ) as 0.5% of R(t, ξ i ). [START_REF] Roques | Using early data to estimate the actual infection fatality ratio from COVID-19 in france[END_REF] conclude that COVID-19 fatalities are under-reported, and using their IFR estimate we adequately fit the trend of the observed COVID-19 deaths but with an offset due to this assumed higher IFR, see Supplementary Materials S1.

Model update Furthermore β estimations can be easily updated as new data become available using parametric empirical Bayes (thus avoiding the need to re-estimate the whole system). It consists in maximizing the likelihood again with respect to β while holding the other parameter distribution fixed to their previously inferred a posteriori distribution. This is how our results are updated with data after March 25 th 2020 in this work.

Effective reproductive number For each region, we compute the effective reproductive number R e (t, ξ i ) as a function of model parameters:

R e (t, ξ i ) = D I b i A(t, ξ i ) + I(t, ξ i ) αA(t, ξ i ) + D q i I(t, ξ i ) D I + D q i (4)
When individuals are homogeneous and mix uniformly, R e (t, ξ i ) is defined as the mean number of infections generated during the infectious period of a single infectious case in the region i at time t. This is the key parameter targeted by NPIs. We compute analytically its 95% confidence interval by accounting for all 95% confidence interval [X min ; X max ] of parameters and trajectories X used in its definition such that:

CI 95% (R e (t, ξ i )) = D I b min i A max (t, ξ i ) + I max (t, ξ i ) αA min (t, ξ i ) + D q min i I min (t, ξ i ) D I + D q max i ; (5) D I b max i A min (t, ξ i ) + I min (t, ξ i ) αA max (t, ξ i ) + D q max i I max (t, ξ i ) D I + D q min i
Asymptomatic proportion At a given time t the number of incident unascertained cases is equal to the sum of two populations, the number of incident non-tested symptomatic individuals (NT) and the number of incident non-tested asymptomatic individuals (AS):

(1 -r)E(t, ξ) D e = N T (t, ξ) + AS(t, ξ) (6)
where r is the proportion of cases tested positive. Collection of data from general practitioners through the re-purposing of the Réseau Sentinelles network to monitor COVID-19 provides a weekly estimation of the number of incident symptomatic cases (tested or not tested) that we previously called r s . This quantity is given over a week but can be evaluated daily by averaging:

r s = rE(t,ξ) De N T (t, ξ) + rE(t,ξ) De ( 7
)
where r s represents the proportion of infected cases seeing a general practitioner. Combining equations ( 6) and ( 7) allows to compute the incident number of asymptomatic cases as a function of the compartment E:

AS(t, ξ) = E(t, ξ) D e 1 -(r + r s ) 1 -r s (8)

Predictions

Short-term predictions of attack rates We predict the proportion of infected individuals among the population in each region at a given date by computing:

[E(t, ξ i ) + A(t, ξ i ) +I(t, ξ i ) + H(t, ξ i ) + R(t, ξ i )] /N i (neglecting the deaths).
Targeting lockdown consequences Given the values of parameters ξ i , we predict the trajectories of the dynamical system compartments using the dlsode differential equation solver in R [START_REF] Soetaert | Solving differential equations in r: Package desolve[END_REF] and we can investigate the impact of NPIs such as lockdown in various scenarios. This impact is driven by two parameters:

• K 2 (K 1 is considered fixed), the decrease ratio of transmission rate of the disease following the first week of NPI, defined as b C i = b i K 2 . This directly translates into a decrease of the effective reproductive number according to Equation (4). It reflects the fact that individual by getting confined decrease their number of contacts. Of note, the current K 2 is estimated by exp(-β 1 -β 2 ) (see Section 2.2.3).

• τ , the duration (in days) of the lockdown during which the mixing and transmission are fixed to b C i = b i K 2 instead of b i . We evaluate the magnitude of the possible epidemic rebound after confinement according to several values for K 2 . In particular, we predict the rates E(t, ξ i )/N i , A(t, ξ i )/N i and I(t, ξ i )/N i on May 11 th 2020 (currently considered by French authorities as the possible start date for lifting lockdown in France). We also compute the optimal lockdown duration τ opt i needed to achieve the epidemic extinction in region i defined as E i (t, ξ i ) < 1 and A i (t, ξ i ) < 1 and I i (t, ξ i ) < 1 simultaneously. In each scenario we predict the date at which the ICU capacities in each region would be overloaded, this date is given by:

t ICU i = arg min ICU (t,ξ i )>η i (t)
η i denoting the ICU capacities limits in region i. We additionally predict how many more ICU beds would be needed at the peak of hospitalization in each scenario, as a proportion of current ICU capacity (DREES, 2019). Finally we provide a rough prediction of the number of deaths for each envisioned scenario assuming a confinement duration of τ opt i days.

3 Results

Estimation of the regional epidemic dynamics

Data fitting Because of the lag inherent to diagnostic testing, we also estimated the number of people already infected at this epidemic start by E 0 (notably, the largest numbers of E 0 in Table 3 are estimated for Île-de-France and Grand Est the two most affected French regions in this early epidemic).

On March 25 th , the cumulative number of ascertained cases was 24,623 and the cumulative number of hospitalized cases was 13,388, see Table 1 for a regional breakdown. Our SEIRAH model fits the data well as can be seen in Figure 2. Moreover, the stability of the estimates is good with a convergence ratio of 1.6 (see Section 2.3.1), corroborating the good identifiability of the estimated parameters (see Supplementary Materials S2). Table 3 provides the regional estimates of the transmission rates (b i ). Of note, regions with higher transmission rate are not necessarily those known to have the highest number of incident ascertained cases. D q i , the number of days from illness onset to hospitalization, can be quite variable between regions and likely accounts for heterogeneity in the observed data. We estimates its population mean at D q 0 = 1.13 days with a standard deviation σ Dq = 0.42.

To evaluate the validity of our structural ODE model ( 1) and of our inference results, we compare the aggregated predictions of the number of both incident ascertained cases and incident hospitalized cases at the national French level to the daily observed incidences (that are still publicly available from SPF at the country level, even after March 25 th -while incident ascertained cases are not openly availilable at the regional level after March 25 th ). Figure 3 displays both predictions and observations, illustrating the added value of incorporating data after March 25 th as those encompass new information about the epidemic dynamics and characteristics as we approach the peak in most regions (notably Grand Est and Île de France). Of note, worse fit of the observed hospitalizations by our model can be explained by the data discrepancy: while we use the SurSaUD R data for our inference (which only account for patients arriving through the emergency room), Figure 3 displays the data from Santé Publique France which should contain all COVID-19 hospitalizations in France (including hospitalizations not coming from the emergency room, as well as Corsica and the French Departements d'Outre Mer that are not taken into account in our model).

q qqqq q qq q qq q q q qq qq q qq q q q qqqqqq qqq q qqqq qq q q q q qq q qqq qq q qq qq qqqq q qq q q q q qq q qqqqq qqqq q q qqq q qq q q q q qq qqqqqqqq q qq qqqq qqq q q q q q q q q q qqqq qqq qqqq qqq qqqq qqqqqq q qq q q q qqqqqqq qqqqqqq qqq q q q qq q qq q q q q q q q q q q q q q qqq qqqq qqq qq q q q qq q q q q qq q qqqq q qq q q q q q q q q q q q q q q q q qqqqqqqq qqqqq qq qqqqq q qq q qqqqqqq q q q q q q q q q q q q q q q q qqq qqq qqq q q q q q qq q q q q q q q qqqqqq q qqqqq qqqq q qq q qq qqqqqqq qqqq q qqq qq q qqq qqqqqqqqqq qq q q q q qqqqqqq qqq q q qqq qqqqqqqq qq qqqq q qq q q qq q q q qqqqqqqq qqq qq qq q qq qq q qqqqqqqq q qq qqqqq q qq q q q q q qqq qqq q q q q q qqq qq q q qq q q q qqqqqqqqqq qqqq q qqqqqqq qqqqqqqqqq qq qqq qqq qq q

Auvergne

Île-de-France Nouvelle-Aquitaine Normandie Occitanie Provence-Alpes-Côte d'Azur Pays de la Loire

Incident confirmed cases

Incident hospitalized cases This estimate is in line with estimation obtained from active surveillance data on the Diamond Princess cruise ship [START_REF] Mizumoto | Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the diamond princess cruise ship, Yokohama, Japan[END_REF].

M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 0 2 M a r 0 9 M a r

Estimation of the lockdown effect

Change of transmission rate during lockdown The parameter β 1 and β 2 measure the effect of the lockdown before and after a week of adjustment. Both are significantly different from 0 (p<0.001) such that the lockdown reduced the transmission rate of COVID-19 by a divisive factor K 1 estimated at 1. 31 [1.27; 1.35] during the first week and K 2 estimated at 3. 63 [3.48; 3.80] after this first week. Of note, thanks to our update algorithm (see Section 2.3.1), it is possible to update those results rapidly as soon as more data are available to inform which scenario of prediction described in Section 3.3 is the most likely.

Effective reproductive number The above quantities directly impact the effective reproductive number as described in

Epidemic dynamics predictions after lockdown lift

In Tables 5 and6, we vary K 2 = 3, 5, 10 (the magnitude of the reduction of transmission during lockdown after the first week). This gradient of simulation is important because the actual French value of K 2 remains currently Table 4: Estimation of the effective reproductive ratios R e during each of the 3 considered periods (before lockdown, during the first week of lockdown, and beyond 1 week of lockdown) for each region with 95% confidence intervals.

unknown. In Section 3.2 we showed that we can estimate a lower bound for K 2 ≥ K 2 = 3.48. From Table 5, we show that the higher K 2 is, the lowest the numbers of ascertained (I), unascertained (A) and latent (E) infected individuals are on May 11 th 2020. However, it is not equal to 0, which means the epidemic is not extinct (and ready to bounce back as soon as lockdown is lifted). In Table 6, we predict the optimal (i.e. shortest) duration of the lockdown to achieve extinction of the epidemic in each region, which is, on average, 407 [236; 786] 5, 97 [88; 112] days if K 2 = 10. We now focus on the average scenario in which K 2 = 5: only three regions in which ICU capacities will be exceeded (Île-de-France, Grand Est and Bourgogne-Franche-Comté) by the end of March 2020. At these times, ICU capacities will have to be increased by respectively 125% [107%; 146%], 178% [147%; 214%] and 158% [125%; 195%]. In this scenario and lockdown is maintained until τ opt i in each region, the predicted number of deaths is 86,387 [71,838; 103,969].

days if K 2 = 3, 147 [123; 179] days if K 2 =
Attack rates under lockdown Table 7 presents the proportion of infected individuals at various dates, i.e. instantaneous attack rates. We also predict them for three horizon dates assuming confinement would be main-tained until these dates: 2020-05-15, 2020-06-08 and 2020-06-22. We predict the national French attack rate on May 15 th 2020 to be 3.8% [3.1%; 4.8%].

Occitanie

Provence-Alpes-Côte d'Azur Pays de la Loire Île-de-France Nouvelle-Aquitaine Normandie

Centre-Val de Loire Grand Est Hauts-de-France 418 [ 3,651; 7,813] 11,479 [ 7,719; 16,720 lockdown lift on May 11 th 2020 We simulated the effect of lifting the lockdown on May 11 th 2020 assuming that after this date the transmission goes back to its value before lockdown. Figure 5 shows the predicted dynamics for each region. In every region, we observed a large rebound occurring either in June or July. The timing and magnitude of this rebound is largely influenced by the importance of the first wave, that is successfully contained thanks to the lockdown. These results strongly argue for enforcing other NPIs when lockdown is lifted in order to contain R e below 1 and prevent this predictable rebound of the epidemic.

Auvergne-Rhône-Alpes Bourgogne-Franche-Comté Bretagne M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 3 0 A p r 0 6 M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 3 0 A p r 0 6 M a r 0 2 M a r 0 9 M a r 1 6 M a r 2 3 M a r 3 0 A p r

Discussion

In this work, we provide estimations of the key parameters of the dynamics of the COVID-19 epidemic in French regions as well as forecasts according to NPIs especially regarding the proportion of infected when lifting the lockdown policy.

The point estimates of the basic reproductive ratios for French regions fluctuated between 2.4 and 3.4 before lockdown took effect, but according to the uncertainty around these estimates they are not substantially different from one region to another. Therefore, observed differences in the number of cases were due to the epicemic starting first in Grand Est and Île-de-France regions. These estimates were close to those reported before isolation using other models [START_REF] Alizon | Report 1: Estimating the basic reproduction number of the covid-19 epidemic in france[END_REF][START_REF] Flaxman | Report 13: Estimating the number of infections and the impact of non-pharmaceutical interventions on covid-19 in 11 european countries[END_REF]. The model provided estimates of the impact of the lockdown on the effective reproductive ratio and although recent data led to a substantial reduction of R e after the lockdown, it remains close to 1 thus without a clear extinction of the epidemic. These estimates should be updated with more recent data that may lead to an estimated R below 1. On the other hand, it is an argument to add other measures such as intensive testing and strict isolation of cases. In addition, the model provides estimates of the size of the population of people who have been or are currently infected. As already reported [START_REF] Domenico | Report #8: Expected impact of school closure and telework to mitigate covid-19 epidemic in france[END_REF], this proportion of subjects is around 2 to 4 percent, so excluding any herd immunity and control of the epidemic by having a large proportion of people already infected and therefore not susceptible. With our estimates of basic reproduction ratio, the epidemic would become extinct by herd immunity with a proportion of 89.5% (95% CI [88.0%; 90.7%]) of infected people.

Interpretation of our results is conditional on the mechanistic model illustrated in Figure 1, and careful attention must be given to the parameters set from the scientific literature and detailed in Table 2, as updated estimates are published every day. First and foremost, our model takes only two kinds of infectious cases into account: confirmed cases I, and unascertained cases A. Our observation model takes I as the number of infectious cases confirmed by a positive PCR SARS-Cov-2 test. Thus, A can be interpreted as unconfirmed symptomatic cases that can be diagnosed by a GP visit (pos-sibly through remote teleconsultation). This is a very simple representation of the COVID-19 infection, which can have various degrees of severity (e.g. asymptomatic, mild, severe) that could be themselves modeled into different compartments. However, very little data is currently available to gather sufficient information to be able to distinguish between those infectious states. Second, our model does not have a compartment for COVID-19 patients in ICU, and the number of occupied ICU beds is simply taken as a fixed percentage (25% based on an estimate from the Bordeaux CHU University Hospital). Meanwhile ICU bed capacity does not account for the recent surge of available ICU beds in response to the COVID-19 epidemic. Compared to [START_REF] Wang | Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan[END_REF], our model does not feature an inflow of susceptibles n (and matching outflow) but population movement across regions are limited during the isolation period (see Supplementary Materials S3 for a thorough discussion). Deaths were also not distinguished from recoveries in the R compartment, but over the observation period this did not impact the main estimates. Third, our model does not take into account the age-structure of the population on the contrary to the recently posted report [START_REF] Salje | Estimating the burden of sars-cov-2 in france[END_REF] using French data. Interestingly, although the models were different and the data not fully identical, our results were comparable. Actually, our approach captures a part of the unexplained variability between regions through the random effects. This variability might be explained at least partly through the difference in age-structure and probability of hospitalization according to the age. We would like to underline the interest of making the data publicly accessible as done by Santé Publique France on the data.gouv.fr web portal, hence allowing our group to work immediately on the topic. Furthermore, we have made our code fully available on GitHub www.github.com/sistm/ SEIRcovid19 and we are currently working on packaging our software for facilitating its dissemination and re-use.

In conclusion, the lockdown has clearly helped controlling the epidemics in France in every region. The number of infected people varies from one region to the other because of the variations in the epidemic start in these regions (both in terms of timing and size). Hence, the predicted proportion of infected people as of May 11 varies, but stays below 10 % everywhere. It is clear from this model, as in other published models [START_REF] Domenico | Expected impact of lockdown in île-de-france and possible exit strategies[END_REF][START_REF] Flaxman | Report 13: Estimating the number of infections and the impact of non-pharmaceutical interventions on covid-19 in 11 european countries[END_REF], that a full and instantaneous lockdown lift would lead to a rebound. Additional measures may help in controlling the number of new infections such as strict case isolation, contact tracing [START_REF] Domenico | Expected impact of lockdown in île-de-france and possible exit strategies[END_REF] and certainly a protective vaccine for which the strategy of administration to the population remains to be defined (Amanat

Figure 1 :

 1 Figure1: SEIRAH model representation -adapted from[START_REF] Wang | Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan[END_REF]

Figure 4 :

 4 Figure 4: Region specific R e (t, ξ i ) compared to the national average. Lockdown started on March 17 th .

Figure 5 :

 5 Figure 5: Region specific predicted dynamics if lockdown is entirely lifted on May 11 th 2020.

Table 1 :

 1 Description of the COVID-19 epidemic by French region. The epidemic start date is defined as the first day followed by three consecutive days with at least 1 confirmed case.

			Observed	Observed			Percentage
			cumulative	cumulative		ICU	of ascertained
		Epidemics Ascertained Hospitalized Population	Capacities	cases from
	Geographical region	start date	cases until	cases until	Size	on 2018-12-31 Réseau Sentinelles
			2020-03-25	2020-03-25			rs
	Auvergne-Rhône-Alpes	2020-03-02	2,082	1,456	8,032,377	559	0.03
	Bourgogne-Franche-Comté	2020-03-02	1,565	1,095	2,783,039	198	0.07
	Bretagne	2020-03-02	602	273	3,340,379	162	0.02
	Centre-Val de Loire	2020-03-10	545	138	2,559,073	180	0.03
	Grand Est	2020-03-02	5,478	2,435	5,511,747	465	0.05
	Hauts-de-France	2020-03-02	1,744	762	5,962,662	438	0.02
	Île-de-France	2020-03-02	7,651	3,870	12,278,210	1,147	0.04
	Nouvelle Aquitaine	2020-03-04	908	695	5,999,982	412	0.03
	Normandie	2020-03-10	677	214	3,303,500	240	0.02
	Occitanie	2020-03-02	1,081	673	5,924,858	474	0.03
	Provence-Alpes-Côte d'Azur	2020-03-02	1,927	1,241	5,055,651	460	0.05
	Pays de la Loire	2020-03-04	363	536	3,801,797	181	0.01
	France	NA	24,623	13,388	64,553,275	4,916	0.03 [0.01;0.07]

Table 2 :

 2 Parameters of the SEIRAH model

Table 3 :

 3 Estimation of region-wise model parameters as well as the final attack rate (without confinement). of the epidemic without intervention It is also interesting to predict the percentage of infected individuals in each region at future dates, corresponding to attack rates. Table3provides such attack rates by region in the absence of any NPI, showing that 89.5% [88%; 90.7%] of the French population would end up being infected. In this scenario we estimate that the peak of the epidemic would have been around May 3 rd 2020 and that the epidemic would have finally gone extinct around May 25 th 2021. At the peak of hospitalization, about714,259 [360,692; 1,065,226] individuals would have required simultaneous hospitalization.Proportion of asymptomatic cases With the current model, we have D

	1 6 M a r 2 3

Figure 2: Fitting curves of I

(in)

and H

(in)

with the SEIRAH model e = 5.1 and r = r s = 3.3% [1.2%; 6.6%]. Thus we estimate that the percentage of asymptomatic infectious individuals is : 18.1% [16.7%; 19.2%].

Table 4

 4

	. The overall ef-

Table 5 :

 5 Epidemic state on May 11 th 2020. For each region depending on K 2 the effect of lockdown on transmission, we present the predicted number of ascertained (I), unascertained (A) and latent (E) individuals infected and its 95% confidence intervals.

			Minimal	Date t ICU i	% of	Number of	Final
		Region	confinement	when ICU	ICU beds	Deaths	attack rates
			duration τ opt i	capacities	occupied	0.005 × R(∞, ξi)	R(∞, ξi)/Ni
			in days	reached	at peak
	K2 = 3	Auvergne-Rhône-Alpes	348 [233; 435]	2020-04-10	107% [ 88%; 127%] 18,567 [15,022; 23,399]	4.6% [3.7%; 5.8%]
		Bourgogne-Franche-Comté	260 [173; 335]	2020-03-25	200% [155%; 254%]	4,045 [ 3,135; 5,388]	2.9% [2.3%; 3.9%]
		Bretagne	228 [156; 286]		62% [ 41%; 88%]	5,080 [ 3,760; 6,891]	3.0% [2.3%; 4.1%]
		Centre-Val de Loire	626 [425; 736]		93% [ 45%; 172%] 20,532 [11,263; 33,550] 16.0% [8.8%; 26.2%]
		Grand Est	366 [256; 447]	2020-03-26	250% [201%; 307%] 26,710 [20,942; 34,837]	9.7% [7.6%; 12.6%]
		Hauts-de-France	259 [184; 314]		61% [ 48%; 76%] 16,211 [13,396; 19,767]	5.4% [4.5%; 6.6%]
		Île-de-France	518 [358; 612]	2020-03-30	205% [164%; 261%] 70,331 [52,721; 98,272] 11.5% [8.6%; 16.0%]
		Nouvelle-Aquitaine	848 [521; 990]	2020-04-13	155% [ 82%; 273%] 32,511 [16,105; 58,929] 10.8% [5.4%; 19.6%]
		Normandie	293 [184; 383]		41% [ 27%; 57%]	7,643 [ 5,546; 10,639]	4.6% [3.4%; 6.4%]
		Occitanie	357 [219; 475]		62% [ 44%; 85%]	9,850 [ 6,861; 14,850]	3.3% [2.3%; 5.0%]
		Provence-Alpes-Côte d'Azur	411 [240; 546]	2020-04-05	129% [ 92%; 172%] 10,524 [ 7,084; 15,864]	4.2% [2.8%; 6.3%]
		Pays de la Loire	373 [232; 478]	2020-04-12	106% [ 66%; 158%] 13,078 [ 8,404; 21,423]	6.9% [4.4%; 11.3%]
	K2 = 5	Auvergne-Rhône-Alpes	148 [119; 170]		79% [ 67%; 93%]	9,185 [ 8,096; 10,529]	2.3% [2.0%; 2.6%]
		Bourgogne-Franche-Comté	124 [ 98; 145]	2020-03-25	158% [124%; 197%]	2,269 [ 1,857; 2,727]	1.6% [1.3%; 2.0%]
		Bretagne	122 [ 97; 142]		52% [ 35%; 73%]	3,229 [ 2,553; 4,063]	1.9% [1.5%; 2.4%]
		Centre-Val de Loire	167 [118; 206]		38% [ 24%; 55%]	3,222 [ 2,405; 4,289]	2.5% [1.9%; 3.4%]
		Grand Est	156 [126; 181]	2020-03-26	178% [148%; 210%] 12,307 [10,584; 14,286]	4.5% [3.8%; 5.2%]
		Hauts-de-France	136 [112; 156]		51% [ 40%; 62%] 10,072 [ 8,755; 11,538]	3.4% [2.9%; 3.9%]
		Île-de-France	181 [144; 209]	2020-03-30	125% [105%; 147%] 23,331 [20,019; 27,317]	3.8% [3.3%; 4.4%]
		Nouvelle-Aquitaine	175 [126; 211]		65% [ 44%; 87%]	4,337 [ 3,172; 5,657]	1.4% [1.1%; 1.9%]
		Normandie	127 [ 98; 151]		31% [ 21%; 42%]	3,933 [ 3,211; 4,818]	2.4% [1.9%; 2.9%]
		Occitanie	141 [110; 166]		45% [ 33%; 61%]	4,553 [ 3,597; 5,782]	1.5% [1.2%; 2.0%]
		Provence-Alpes-Côte d'Azur	146 [113; 171]		88% [ 68%; 112%]	4,232 [ 3,423; 5,231]	1.7% [1.4%; 2.1%]
		Pays de la Loire	146 [111; 172]		74% [ 49%; 105%]	5,717 [ 4,166; 7,732]	3.0% [2.2%; 4.1%]
	K2 = 10 Auvergne-Rhône-Alpes	102 [ 87; 114]		72% [ 60%; 84%]	7,219 [ 6,431; 8,034]	1.8% [1.6%; 2.0%]
		Bourgogne-Franche-Comté	89 [ 75; 101]	2020-03-25	144% [115%; 176%]	1,838 [ 1,547; 2,167]	1.3% [1.1%; 1.6%]
		Bretagne	90 [ 75; 102]		48% [ 32%; 67%]	2,701 [ 2,147; 3,292]	1.6% [1.3%; 2.0%]
		Centre-Val de Loire	93 [ 74; 109]		30% [ 19%; 45%]	2,079 [ 1,631; 2,615]	1.6% [1.3%; 2.0%]
		Grand Est	105 [ 90; 118]	2020-03-26	157% [131%; 186%]	9,426 [ 8,149; 10,859]	3.4% [3.0%; 3.9%]
		Hauts-de-France	99 [ 85; 111]		47% [ 38%; 57%]	8,353 [ 7,392; 9,444]	2.8% [2.5%; 3.2%]
		Île-de-France	114 [ 98; 128]	2020-03-31	107% [ 92%; 124%] 16,766 [14,859; 19,112]	2.7% [2.4%; 3.1%]
		Nouvelle-Aquitaine	101 [ 82; 116]		53% [ 37%; 72%]	2,856 [ 2,168; 3,665]	1.0% [0.7%; 1.2%]
		Normandie	87 [ 71; 100]		28% [ 19%; 38%]	3,104 [ 2,588; 3,721]	1.9% [1.6%; 2.3%]
		Occitanie	97 [ 80; 110]		41% [ 30%; 54%]	3,538 [ 2,825; 4,463]	1.2% [1.0%; 1.5%]
		Provence-Alpes-Côte d'Azur	97 [ 81; 111]		78% [ 59%; 98%]	3,202 [ 2,576; 3,901]	1.3% [1.0%; 1.5%]
		Pays de la Loire	98 [ 79; 111]		66% [ 43%; 92%]	4,361 [ 3,158; 5,763]	2.3% [1.7%; 3.0%]

Table 6 :

 6 Predicted outcomes under lockdown during optimal time τ opt i . Note that ICU and Deaths are not directly modeled as a compartment in our model and based on crude estimates. ICU capacities are considered as of 2018-12-31 according to DREES (2019) and do not account for the recent surge in ICU beds in France.

Acknowledgements

The authors thank Romain Griffier for his time in discussing the aggregated features of COVID-19 patients care at the Bordeaux University Hospital. The authors thank the opencovid-19 initiative for their contribution in opening the data used in this article. BPH thanks Vincent Pey for discussions about the clinical characteristics of the COVID-19 infection. This work is supported in part by Inria Mission COVID19, project GESTEPID.

Availability

The data from the SurSaUD R database regarding COVID-19 is available from the data.gouv French government platform at https://www.data. gouv.fr/fr/datasets/donnees-des-urgences-hospitalieres-et-de-sosmedecins-relatives-a-lepidemie-de-covid-19. The source code used for this work is available on GitHub at www.github.com/sistm/SEIRcovid19.

 [START_REF] Amanat | Sars-cov-2 vaccines: status report[END_REF][START_REF] Lurie | Developing Covid-19 Vaccines at Pandemic Speed[END_REF][START_REF] Thanh | The COVID-19 vaccine development landscape[END_REF].

Competing interest

The authors have no competing interests to declare.

Author contributions

MP, LW, RT and BPH designed the study. MP and BPH analyzed the data. MP, DD and BPH implemented the software code. QC performed identifiability and asymptotic analysis of the model. MP, LW, RT and BPH interpreted the results and wrote the manuscript.