M U B Niazi

C Canudas-De-Wit
email: carlos.canudas-de-wit@gipsa-lab.fr

A Kibangou
email: alain.kibangou@univ-grenoble-alpes.fr

Optimal Control of Urban Human Mobility for Epidemic Mitigation Muhammad Umar B. Niazi

Carlos Canudas-de-Wit Alain Kibangou Pierre-Alexandre Bliman

Abstract-Controlling human mobility during an epidemic is a fundamental issue faced by policymakers. Such control can only be done optimally if human mobility is adequately modeled at the scale of a city or metropolis. This paper, first, develops a model of human mobility that captures the daily patterns of mobility in an urban environment through time-dependent gating functions, which are controlled by the destination schedules and mobility windows. The process of epidemic spread is incorporated at each location that depends on the number of susceptible and infected people present at that location. Then, two optimal control policies are proposed to maximize the economic activity at the destinations while mitigating the epidemic. Precisely, operating capacities and time schedules of destinations are controlled to maximize the economic activity under the constraint that the number of active infected cases remains bounded.

I. INTRODUCTION

Although human mobility plays a vital role in a country's economy, it facilitates the spread of disease in the event of an epidemic by allowing contact between infected and susceptible populations. If human mobility is not controlled during an epidemic, it can result in a huge number of infected cases, which could overwhelm the hospitals and cause the loss of lives. On the other hand, strict restrictions on human mobility can also halt the economy and result in the loss of livelihoods. Therefore, the models of human mobility that incorporate the process of disease spread play a vital role in the analysis, understanding, and mitigation of epidemics.

Considering a SIR epidemic model for disease spread, human mobility between different geographic regions has been investigated and modeled in [START_REF] Sattenspiel | A structured epidemic model incorporating geographic mobility among regions[END_REF]- [START_REF] Poletto | Human mobility and time spent at destination: Impact on spatial epidemic spreading[END_REF]. In these models, the individuals associated with one region or city can go to and return from other regions. However, these models capture averaged mobility patterns between different cities with large timescales and cannot capture the daily patterns of mobility within an urban environment. To tackle this problem, [START_REF] Frias-Martinez | An agentbased model of epidemic spread using human mobility and social network information[END_REF]- [START_REF] Nadini | A multi-agent model to study epidemic spreading and vaccination strategies in an urban-like environment[END_REF] study agent-based models of urban human mobility with epidemic spread. These models are powerful tools for computational purposes, however, they rely on the digital footprints of individuals and can lead to privacy violations. Another line of research in [START_REF] Song | Reinforced epidemic control: Saving both lives and economy[END_REF] relates to the control of human mobility. It aims to find control policies that restrict mobility to and from regions that are estimated to be of high risk by employing reinforcement learning framework and relying on aggregated demand for mobility and regional epidemic statistics.

In this paper, we develop an urban human mobility model that captures the daily mobility patterns and incorporates the process of epidemic spread at each location. Every day a certain number of people go from their residential areas, which are called origins, to locations visited daily for work, education, shopping, etc., which are called destinations, and return on the same day. The daily mobility patterns are captured by the time-dependent supply and demand gating functions. The supply gating function (SGF) of each destination is controlled by its daily destination schedule, which is its opening and closing hours. The demand gating function (DGF), on the other hand, is defined on each edge of the mobility network and corresponds to the daily mobility window, which is the time interval during which people utilize that edge to move between origins and destinations. The supply function of each destination controlled by the SGF determines the inflow allowed to that destination and depends on its operating capacity controlled by the capacity control input. The demand function controlled by the DGF determines the outflow from one location to another. The process of urban human mobility is modeled on the network edges that connect different locations through flows and the process of epidemic spread is modeled locally at each location that depends on the number of susceptible and infected people at that location.

We formulate two optimal control problems for epidemic mitigation while maximizing the economic activity: (i) optimal capacity control policy and (ii) optimal schedule control policy. These problems aim to find an optimal capacity control input and schedule control input, respectively, that maximizes the economic activity while mitigating the epidemic by keeping the number of active infected cases bounded. The capacity control policy restrict the number of people in destinations of each category by specifying the operating capacities in relation to their nominal capacities. The schedule control inputs, on the other hand, specifies the closing hour of destinations of each category by altering the destination schedules and mobility windows.

The paper is organized as follows: Section II introduces the urban human mobility model that incorporates the process of epidemic spread at each location. Section III presents optimal capacity control and optimal schedule control policies, and illustrates their efficacy through a numerical example. Finally, Section IV presents the concluding remarks.

II. MODEL OF URBAN HUMAN MOBILITY WITH EPIDEMIC SPREAD

Consider human mobility in an urban environment between locations of two types: origins and destinations. The origins correspond to locations where people reside-for example, residential areas, neighborhoods, and towns. The destinations, on the other hand, correspond to locations that people visit daily for work, education, shopping, or leisurefor example, industrial zones, business parks, schools, markets, cinemas, etc. Every day a certain number of people go from each origin to the destinations during specified time intervals and then return later the same day. At each location, there is an epidemic spread process that will be described here by a SIR model, [START_REF] Hethcote | The mathematics of infectious diseases[END_REF], where the population is divided into Susceptible, Infected, and Recovered classes, and the disease is transmitted according to the local infection rates when the susceptible and infected populations mix in the same location. Notice that any similar epidemiological model could be used instead.

A. Network Representation and Main Assumptions

Let V o = {1, . . . , m} be the index set of m origins and V d = {m + 1, . . . , m + n} be the index set of n destinations. The network of urban human mobility is represented by a bi-directed, bipartite graph G = (V o , V d , E), where E is the set of bi-directed edges-i.e., for every i ∈ V o and j ∈ V d , if (i, j) ∈ E then (j, i) ∈ E. Figure 1 illustrates a mobility network, which is the running example of this paper.

We adopt the following assumptions in the model: (A1) The total population of the city remains constant. (A2) The mobility occurs only between pairs of origins and destinations, and not among a pair of different origins or a pair of different destinations. (A3) The number of people who visit destination j from origin i during a day is equal to the number of people who return to i from j on the same day. (A4) The mobility pattern between each pair of origins and destinations is periodic and repeats every day, i.e., the period T period = 24 hours. In particular, we ignore mobility patterns of the weekends or holidays that may be different than the normal days. (A5) The restrictions imposed on the urban human mobility by the government affects all the people, whether susceptible, infected, or recovered, equally. Note that (A4) is a simplifying assumption, which can be relaxed without loss of generality by considering T period to be a week, a month, etc.

Denote P i to be the total population of origin i ∈ V o , which is the number of people who reside in i, and C j to be the nominal capacity of destination j ∈ V d , which is the maximum number of people who can visit j at one time when there is no epidemic. Note that the total population of the city P = i∈Vo P i .

B. Destination Categories

Suppose the destinations are divided into p ≤ n categories, which correspond to workplace, school, market, etc. The In the example illustrated in Figure 1, we have p = n, i.e., three categories and three destinations, because each destination is of a different category. However, it is possible that a mobility network may have multiple destinations of the same category. For instance, a slight modification of the example of Figure 1, say with two working places and three schools, could illustrate this immediately.

C. Operating Capacities of Destinations

Let u k (t) ∈ [0, 1] be the capacity control input for destinations of category k, for k = 1, . . . , p, which determines the allowed operating capacity of D k in terms of the proportion of nominal capacity at time t in the event of an epidemic. In other words, it can be considered as the government's policy at time t that limits the operating capacities in destinations of category k in order to mitigate the epidemic spread, where

Operating capacity = C j u k (t), for j ∈ D k .
We consider u k (t) to be piece-wise constant-i.e.,

u k (t) =          µ 1 k if t ∈ [0, T u ) µ 2 k if t ∈ [T u , 2T u ) . . . . . . µ q k if t ∈ [(q -1)T u , qT u ) (1)
with µ h k ∈ [0, 1] constant for every h ∈ {1, . . . , q} and T = qT u the total time horizon considered by a policymaker. The policy horizon T u is a multiple of T period and corresponds to the time after which the policy on operating capacities is announced periodically. It can be chosen by the policymaker from at least a week to several months because changing the policy on shorter time intervals may not be practical in terms of implementation.

D. Destination Schedule and Mobility Window

The destination schedule of j ∈ V d

S j = [a j , b j ), 0 ≤ a j < b j ≤ 24
is the daily time interval during which j is open, where a j and b j are the nominal opening and closing hours of j, respectively. The origins, on the other hand, are open throughout the day, i.e., for every i ∈ V o , S i = [0, 24).

Then, the supply gating function

(SGF) for j ∈ V o ∪ V d σ j (t) = 1 if t mod 24 ∈ S j 0 otherwise (2)
which is periodic with respect to 24 hours. The mobility window of (i, j) ∈ E

D ij = [t ij , t ij + τ ij ), 0 ≤ t ij < t ij + τ ij ≤ 24
is the daily time interval during which there is mobility from i to j, where τ ij > 0 is the duration of mobility window in hours. Then, the demand gating function (DGF) of (i, j) ∈ E

δ ij (t) = 1 if t mod 24 ∈ D ij 0 otherwise (3)
which is also periodic with respect to 24 hours.

E. Model of Urban Human Mobility

Let N i (t) ≥ 0 be the number of people in i ∈ V o ∪ V d at time t (hour). Then, according to the urban human mobility model, the rate of change of the number of people at any location at time t is equal to the sum of inflows to that location minus the sum of outflows from that location.

In other words, for any i ∈ V o ∪ D k and k ∈ {1, . . . , p}, the urban human mobility model is given by

Ṅi = j∈Ni (φ ji -φ ij ) (4)
where N i is the set of neighbors of i in the mobility network

G and φ ij (t, N i (t), N j (t), u k (t)) is the flow from i ∈ V o to j ∈ D k given as φ ij = min(∆ ij , Σ j ) with ∆ ij (t, N i (t), u k (t)
) and Σ j (t, N j (t), u k (t)) the demand and supply functions, respectively. Notice that the flow φ ji (t, N j (t), N i (t), u k (t)) is defined similarly with the subscript ji instead of ij.

The supply function Σ j (t, N j (t), u k (t)) of each location j corresponds to the allowed inflow to j from other locations and is given by

Σ j = σ j min(F j , v[C j u k -N j ]) if j ∈ D k σ j min(F j , v[P j -N j ]) if j ∈ V o i j φ ij (t) = min(∆ ij , Ψ j ) Σ j = σ j min(F j , v[C j u k -N j ]) Ψj ∆ ij = δ ij min(vN i , f ij ) ∆ij P i f ij N i (t) ∆ ij Slope v C j F j N j (t) Ψ j Slope -v
Fig. 2: An example illustrating the flow φ ij (t) from origin i to destination j in terms of supply of j and demand of i with respect to j. Here, the arrows on each curve indicate the time evolution.

where σ j (t) is the SGF given by (2), v > 0 is a regularization parameter taken to be very large (see Remark 1), C j u k (t) is the operating capacity of j ∈ D k with u k (t) defined in (1), and

F j (t) = i∈Nj f ij (t)
is the maximum inflow to j with

f ij (t) = M ij u k (t) τ ij (5)
the maximum outflow from i ∈ V o to j ∈ D k . Here, M ij denotes the nominal number of visitors to j from i and M ij u k (t) is the number of visitors when the capacity control input u k (t) is implemented. Notice that, since all visitors return, we have

M ji = M ij . The demand function ∆ ij (t, N i (t), u k (t))
of each edge (i, j) ∈ E corresponds to the outflow from i towards j and is given by

∆ ij = δ ij min (vN i , f ij )
where δ ij (t) is the DGF given by (3), v > 0 is the same regularization parameter introduced in the supply function Σ j (t, N j (t), u k (t)), and f ij (t) is the maximum outflow from i to j given by [START_REF] Frias-Martinez | An agentbased model of epidemic spread using human mobility and social network information[END_REF]. Suppose i ∈ V o and j ∈ D k , then Figure 2 illustrates the flow from i to j in terms of demand and supply functions. In the figure, notice that the demand of i moves from right to left with respect to time t, i.e., from being full to being empty, and the supply of j moves from left to right with respect to time t, i.e., from being empty to being full, which is indicated by arrows in the figure .   Remark 1. To ensure that the daily number of people going from i to j equals M ij u k (t), we assume that the demand function

∆ ij (t, N i (t), u k (t)) = δ ij (t)f ij (t)1 Ni(t)>0
, where 1 Ni(t)>0 = 1 if N i (t) > 0, and 0 otherwise, is the indicator function. Similarly, we assume that the supply function Ψ j (t, N j (t), u k (t)) = σ j (t)F j (t)1 Nj (t)<Cj u k (t) . However, to avoid the discontinuity posed by the indicator functions, we approximate the demand and supply functions by considering steep slope with a very large regularization parameter v as illustrated in Figure 2.

Capacity control policy u 1 (t) = 0.5, u 2 (t) = 0.5, u 3 (t) = 1 for t ∈ [0, 48) Regularization parameter v = 30
TABLE I: Parameters related to urban human mobility for the example of Figure 1.

Example 1. Consider the example of mobility network shown in Figure 1 with two origins and three destinations. For the mobility model (4), we consider the parameters given in Table I. The mobility profile of two days is plotted in Figure 3, where N 1 (t), N 2 (t) denote the number of people in Town 1 and 2 at time t, respectively, and N 3 (t), N 4 (t), N 5 (t) denote the number of people in the industries, schools, and markets, respectively. As shown in the figure, people go from the origins (1 and 2) to the destinations (3, 4, and 5) and return on the same day according to the destination schedules and mobility windows. Notice that the mobility profiles are the same for both days because, by assumption (A4), the destination schedules and mobility windows are the same for every day. Moreover, the nominal capacities of industries and schools are C 3 = 2000 and C 4 = 1500, however, the capacity control policy u 1 (t) = u 2 (t) = 0.5 reduce the operating capacities to 50% of the nominal capacities. Therefore, the maximum number of people present in these destinations during a day is around C 3 /2 = 1000 and C 4 /2 = 750, respectively.

F. Incorporating the Epidemic Spread Process

We consider a susceptible, infected, recovered (SIR) model of epidemic spread, which divides the number of people N i (t) at each location i ∈ V o ∪ V d into three classes: number of susceptible S i (t), infected I i (t), and recovered R i (t), where, at every time t, Fig. 3: The mobility profile of two days for the example in Figure 1.

N i (t) = S i (t) + I i (t) + R i (t).
The disease transmission at each location i occurs according to the local mass action law

β i (t)S i (t) I i (t) N i (t)
where

β i (t) = β i Ni(t) Pi if i ∈ V o β i Ni(t) Ci if i ∈ V d
is the infection rate of i at time t with β j the nominal infection rate of i. The nominal infection rate is defined as the average number of contacts of a person in location i per hour when the number of people in i is maximum. The infection rate β i (t) reduces when the number of people N i (t) at location i is small and increases when N i (t) is large. The infected people I i (t) recover with a recovery rate γ ∈ (0, 1], which is a constant that depends on the disease biology and, if available, the treatment methods. The recovery rate γ is defined as the inverse of the average recovery period (in hours) of the infected cases.

Recall the assumption (A5), which implies that all people irrespective of whether they are susceptible, infected, or recovered can go from one location to another with equal probability. That is, the flow from i to j in terms of the number of susceptible, infected, and recovered can be respectively given by

φ ij (t, N i (t), N j (t), u k (t)) Si(t) Ni(t) φ ij (t, N i (t), N j (t), u k (t)) Ii(t) Ni(t) φ ij (t, N i (t), N j (t), u k (t)) Ri(t) Ni(t) . Let x i (t) = [ S i (t) I i (t) R i (t) ] T ∈ R 3 ≥0 be the state vector of location i ∈ V o ∪ V d and ξ i (x i (t)) =    -β i (t)S i (t) Ii(t) Ni(t) β i (t)S i (t) Ii(t) Ni(t) -γI i (t) γI i (t)    ∈ R 3 (6) S i I i R i βi γ S j I j R j βj γ φ ij φ ji Origin i Destination j
Fig. 4: The process of urban human mobility happens along the edges whereas the process of epidemic spread happens inside the locations.

be the vector describing the process of epidemic spread in location i. Then, for i ∈ V o and j ∈ D k , the model of urban human mobility with epidemic spread is given by

ẋi = ξ i (x i ) + p k=1 j∈D k φ ji (u k ) x j N j -φ ij (u k ) x i N i ẋj = ξ j (x j ) + i∈Vo φ ij (u k ) x i N i -φ ji (u k ) x j N j . (7)
As illustrated in Figure 4, there are two aspects of the model. First, inside the locations i and j, there is a process of epidemic spread that transmits the disease from the infected to the susceptible with a local infection rate, and the recovery process that heals the infected with a constant recovery rate. Second, on the edges (i, j) and (j, i), there is a process of human mobility that transfers people from one location to another through the flows φ ij (t, N i (t), N j (t), u k (t)) and φ ji (t, N j (t), N i (t), u k (t)), respectively.

III. OPTIMAL MOBILITY CONTROL FOR EPIDEMIC MITIGATION

A. Compact Representation of the Model

The urban human mobility model with epidemic spread given in [START_REF] Pastor-Satorras | Epidemic processes in complex networks[END_REF] describes the dynamics of the number of susceptible S i (t), infected I i (t), and recovered R i (t) people in location i ∈ V o ∪V d . These dynamics are controlled by the piece-wise constant capacity control inputs u 1 (t), . . . , u p (t), defined in (1), of the p destination categories in D. Define

u(t) = u 1 (t) . . . u p (t) T .
Then, to represent the model in a compact form, let

x = x o x d
, where

x o = [ x T 1 • • • x T m ] T x d = [ x T m+1 • • • x T m+n ] T with x o (t) ∈ R 3m ≥0
the state vector of origins V o and x d (t) ∈ R 3n

≥0 the state vector of destinations V d . Similarly, let

ξ = ξ o ξ d , where ξ o = [ ξ T 1 • • • ξ T m ] T ξ d = [ ξ T m+1 • • • ξ T m+n ] T
with ξ o (x o (t)) ∈ R 3m the vector describing the epidemic process in the origins V o and ξ d (x d (t)) ∈ R 3n the vector describing the epidemic process in the destinations V d . Notice that the vector ξ i (x i (t)) describing the epidemic process in each location i is given in [START_REF] Pappalardo | Returners and explorers dichotomy in human mobility[END_REF].

The model ( 7) can be represented as

ẋ = ξ(x) + Φ(x, u)x (8)
where the dependence on t is omitted for brevity. The matrix of flows Φ(t, x(t), u(t)) describes the mobility process in the network G and is given as

Φ(x, u) = Φ oo (x, u) Φ do (x, u) Φ od (x, u) Φ dd (x, u) ⊗ I 3
where ⊗ denotes the Kronecker product and

Φ oo = diag - j∈V d φ 1j N 1 . . . - j∈V d φ mj N m Φ dd = diag - i∈Vo φ m+1,i N m+1 . . . - i∈Vo φ m+n,i N m+n Φ od =       φ 1,m+1 N 1 • • • φ m,m+1 N m . . . . . . . . . φ 1,m+n N 1 • • • φ m,m+n N m       Φ do =       φ m+1,1 N m+1 • • • φ m+n,1 N m+n . . . . . . . . . φ m+1,m N m+1 • • • φ m+n,m N m+n       .

B. Economic Activity and Active Infected Cases

The economic activity E(t) ∈ R ≥0 in the mobility network G at time t is defined as

E(t) = p k=1 j∈D k χ j N j (t) C j (9)
where N j (t) is the number of people in destination j ∈ D k at time t, C j is the nominal capacity of j, and χ j ∈ [0, 1] is the weight assigned to j according to its economic importance such that j∈V d χ j = 1. Since N i (t) = 1 T 3 x i (t), we can write [START_REF] Song | Reinforced epidemic control: Saving both lives and economy[END_REF] as

E(t) = e T x(t) (10)
where x(t) is the state of ( 8) at time t and

e = 0 T m χm+1 Cm+1 • • • χm+n Cm+n T ⊗ [ 1 1 1 ] T .
The number of active infected cases I(t) ∈ R ≥0 in the mobility network G at time t is the sum of the number of infected people in all the locations. It is given by

I(t) = i∈Vo∪V d I i (t) (11)
where I i (t) is the number of infected people at location i at time t. Since I i (t) = [ 0 1 0 ]x i (t), we can write (11) as where

I(t) = g T x(t)
g = 1 T m 1 T n T ⊗ [ 0 1 0 ] T .
Finally, the infection peak is defined as

I peak = sup t∈[0,T ] I(t) = sup t∈[0,T ] g T x(t) (12)
where [0, T ] is a given finite time horizon. Example 2. Consider the example of network in Figure 1 with the mobility parameters given in Table I and the epidemic parameters given in Table II. Note that the nominal infection rates outside the residences (the "destinations") are assumed higher than the nominal infection rates at the residences (the "origins"). For T = 1680 hours (or 10 weeks), we plot the active infected cases I(t) in Figure 5 under two circumstances: (i) when there are no restrictions on the operating capacities of all destinations, i.e., u(t) = [ 1 1 1 ] T for all t ∈ [0, T ], and (ii) when the operating capacities of all destinations are reduced to 50% throughout [0, T ], i.e., u(t) = [ 0.5 0.5 0.5 ] T for all t ∈ [0, T ]. In Figure 5, notice the effect on the number of active infected cases when the operating capacities are reduced. When there are no restrictions on the operating capacities, the infection peak I peak is about 1810 people, whereas, with the restrictions, the peak is about 605 people.

In other words, given the mobility and epidemic parameters for the example of Figure 1, one can reduce the infection peak by 66.5% through 50% reduction of the operating capacities. TABLE III: Parameters related to the optimal control problem (13) for the example of Figure 1.

C. Optimal Capacity Control Policy

The number of hospitalized cases and deaths due to an epidemic are related to the number of infected cases I(t). A large value of I(t) implies a large number of hospitalizations and loss of lives in the near future. Moreover, if the measures to mitigate the epidemic and limit the number of infected cases are not taken, then the number of hospitalizations may reach a point that could challenge the available medical facilities of the city. As shown in Example 2, the number of infected cases can be reduced by reducing the operating capacities of the destinations through the capacity control policies u 1 (t), . . . , u p (t). However, on the other hand, choosing the values of u 1 (t), . . . , u p (t) too small can result in a significant reduction of the economic activity E(t), which may result in bankruptcy of businesses and loss of livelihoods. Therefore, our goal in this section is to find optimal capacity control policy that maximize the economic activity under a constraint that the infection peak I peak remains bounded from above.

Suppose a finite time horizon T , a policy horizon T u = T /q for q ∈ N, an upper bound I > 0 on the infection peak, and all the parameters of the model (8) be given. Let u ∈ U, where U is the set of admissible capacity control policies u : [0, T ] → [0, 1] p such that, for every t ∈ [(h -1)T u , hT u ) and h ∈ {1, . . . , q}, u(t) = µ h for some

µ h = [ µ h 1 • • • µ h p ] T ∈ [0, 1] p .
Then, the optimal capacity control policy is obtained by solving the following problem:

maximize u∈U L(u) := 1 T T 0 e T x(t; x 0 , u)dt subject to ẋ = ξ(x) + Φ(x, u)x; x(0) = x 0 I peak (u) ≤ I (13)
where I peak (u) = sup t∈[0,T ] g T x(t; x 0 , u) is given in ( 12) and the economic activity E(t) = e T x(t) is given in [START_REF] Hethcote | The mathematics of infectious diseases[END_REF]. Example 3. Again, consider the example of Figure 1 with the mobility and epidemic parameters given in Table I and II, respectively. Also, consider the parameters in Table III required by the optimal control problem (13). The optimal control problem is solved numerically for a time horizon T = 1680 hours (i.e., 10 weeks) using a nonlinear programming solver fmincon in MATLAB with interior point algorithm. The solver returns a local minimum u opt (t) plotted in Figure 6 that satisfies the constraints (??) and is piecewise constant, where the policy horizon T u = 336 hours (i.e., 2 weeks). In particular, the constraint on the infection peak is satisfied and I peak ≈ 922 is less than I = 1000 as shown in Figure 7. Notice that in the beginning the optimal capacity control allows the operating capacities to be around 70-90% of the nominal capacities of destinations. However, as the number of infected cases increase, the value of the optimal capacity control decreases for the next four steps until 8T u = 1344 hours (i.e., 8 weeks) to mitigate the infection spread. Then, in the last interval [8T u , T ], the optimal capacity control increases to allow more people visiting the destinations because the infected cases have started to decrease. In Figure 8, we plot the economic activity for three cases: (i) when no capacity control policy is implemented, (ii) when the capacity control policy limits the operating capacities to 50% of the nominal capacities of destinations, and (iii) when the optimal capacity control policy u opt (t) shown in Figure 6 is applied. Notice that the optimal capacity control u opt (t) increases the economic activity as compared to the cases when u(t) = 0.51 3 while keeping the I peak under the bound I. Fig. 8: Economic activity (i) without capacity control policy, i.e., u = 1 3 , (ii) with 50% capacity control policy u = 0.51 3 , and (iii) with optimal capacity control policy u opt .

D. Optimal Schedule Control Policy

We formulate an epidemic mitigation policy that alters the destination schedules and mobility windows. That is, for every j ∈ D k , the destination schedule is altered as

S j = [a j , min(s k , b j ))
where s k ∈ [s, 24) is the schedule control that enforces that all destinations of category k, for k = 1, . . . , p, must be closed after s k hour, respectively, and s ≥ 0 is the lower bound on s k . Such a policy limits the spread of infection by reducing the daily amount of time people spend at destinations. It also alters the mobility windows

D ij = [t ij , min(t ij + τ ij , s k )) D ji = [min(t ji , s k ), min(t ji , s k ) + τ ji ).
for i ∈ V o and j ∈ D k . That is, people cannot go from i to j after s k hour and people at j must return to i after s k hour.

Let s = [ s 1 • • • s p ] T be the schedule control policy of all destination categories. Then, the problem is to find an optimal s such that the economic activity is maximized while keeping the infection peak bounded by I. The schedule control policy complements the capacity control policy obtained by solving (13) when there are lower bounds on the capacity control policy. These lower bounds correspond to minimum operating capacities of certain destinations that are required for functioning of the society. This is because some destinations, like hospitals and markets, are essential and their operating capacities cannot be reduced beyond a minimum bound. In other words, for all k ∈ {1, . . . , p} and h ∈ {1, . . . , q}, we assume constant capacity control policy u(t) = µ, for all t ∈ [0, T ], where µ ∈ [0, 1] p states the minimum allowed capacity control policy of p destination categories. In the presence of these lower bounds, the problem (13) may become infeasible and the infection peak may no longer be bounded. Thus, implementation of optimal schedule control policy s may help in containing the infections while also allowing economic activity at destinations. Suppose a finite time horizon T , an upper bound I > 0 on the infection peak, a constant capacity control policy u(t) = µ ∈ [0, 1] p , ∀t ∈ [0, T ], and the parameters of the model (8) are given. Then, the optimal schedule control policy is obtained by solving the following problem: where I peak (s) = sup t∈[0,T ] g T x(t; x 0 , s) is given in (12) and the economic activity E(t) = e T x(t) given in [START_REF] Hethcote | The mathematics of infectious diseases[END_REF].

Example 4. Consider again the example of Figure 1 and the parameters given in Table I, II, and III, where we change the policy horizon T u = T since the schedule control policy s is constant throughout [0, T ]. Let the capacity control policy u(t) = [ 0.75 0.75 0.75 ] T , ∀t ∈ [0, T ], i.e., the operating capacities of destinations should be 75% of their nominal capacities. Then, solving (14) numerically by using fmincon solver in MATLAB with interior point algorithm, we obtain a local minimum s opt = [ 14.54 13.75 16.10 ] T , which means that the industries must be closed after 14.54 hour (02:32 pm), the schools must be closed after 13.75 hour (01:45 pm), and the markets must be closed after 16.10 hour (04:06 pm). Figure 9 shows that the constraint on the infection peak is satisfied.

IV. CONCLUSION

We developed an urban human mobility model on network of origins and destinations that incorporates the process of epidemic spread at each location. The model is described by the flows that transfer people from origins to destinations and back to origins every day. The flows capture the daily patterns of mobility in an urban environment through the gating functions that depend on the destination schedules and mobility windows. At each location, the disease spreads through the interaction of susceptible and infected people, where the infection rate depends on the number of people in that location. We study two optimal control policies, capacity control and schedule control, that mitigate the epidemic spread while maximizing the economic activity at each destination. The optimal capacity control policy maximizes the economic activity by allowing maximum allowable number of people at each destination under the constraint that the infection peak remains bounded. The optimal schedule control policy maximizes the economic activity by allowing maximum allowable time that people spend daily at destinations under the constraint that the infection peak remains bounded.

The future investigations include the model predictive control formulation of optimal control problems. A work in progress is the validation of model with the urban mobility data and regional epidemic statistics and the development of web interface to illustrate our proof of concept.

Fig. 1 :

 1 Fig.1: An example of an urban human mobility network with two origins and three destinations.

Populations P 1

 1 = 3000, P 2 = 2000 Capacities C 3 = 2000, C 4 = 1500, C 5 = 200 No. of Visitors M 13 = 1200, M 14 = 900, M 15 = 900 M 23 = 800, M 24 = 600, M 25 = 600 Mobility windows D i3 = [8, 9.5), D i4 = [9, 10), D i5 = [8.5, 20) D 3i = [17, 18.5), D 4i = [16, 17), D 5i = [10, 21) for i = 1, 2 Destination schedules S 3 = [8, 18), S 4 = [9, 17), S 5 = [8.5, 20)

Recovery rate γ = 1 /Fig. 5 :

 15 Fig. 5: Reduction of the number of infected cases I(t) by reducing operating capacities via capacity control policy.

χ 3

 3 = 0.4, χ 4 = 0.3, χ 5 = 0.3

Fig. 6 :Fig. 7 :

 67 Fig. 6: Optimal capacity control policy u opt .

Fig. 9 :

 9 Fig.9: Reduction of the number of active infected cases I(t) by controlling destination schedules and mobility windows via optimal schedule control policy s opt .

e

  T x(t; x 0 , s)dt subject to ẋ = ξ(x) + Φ(x, s)x; x(0) = x 0 I peak (s) ≤ I (14)

This work is partially supported by European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERCAdG no. 694209, Scale-FreeBack, website: http:// scale-freeback.eu/) and by Inria, France, in the framework of the Inria Mission Covid-19 (Project HealthyMobility).


associate professor in Agnès Caille
email: agnes.caille@med.univ-tours.fr

reader in Sally Kerry

PhD Elsa Tavernier

Clémence Leyrat

professor of biostatistics Sandra Eldridge

Bruno Giraudeau

Timeline cluster: a graphical tool to identify risk of bias in cluster randomised trials

 Cluster randomised trials can be at risk of bias when participants are identified and recruited after randomisation  Reports of cluster randomised trials often fail to adequately describe the recruitment process and whether participants and trial staff are blinded to allocation status at key stages of the trial  This article presents a graphical tool depicting the time sequence and blinding status of the different stages of a cluster randomised trial together with examples to help researchers describe the storyline of such trials  Our graphical tool should be used at both the protocol and report-writing stages to clarify the trial process and to help identify the risk of bias

Standfirst

Robust evidence of the effectiveness of interventions relating to policy, practice and organisation of healthcare often comes from well-conducted cluster randomised trials.

However, such trials are prone to recruitment bias depending on whether participants are recruited before the randomisation of clusters and whether the recruiter is blinded to the allocation status. In most cases, participants and trial staff cannot be blinded to the intervention, which may lead to performance and detection bias. Unfortunately, cluster trial reports often do not provide a clear description of the timing of trial processes and blinding, and these aspects are not covered by current reporting tools. This article proposes a graphical tool depicting the time sequence of steps and blinding status in cluster randomised trials. The tool may be helpful at both the protocol and report-writing stages to clarify the process and to help identify potential bias in cluster randomised trials.

31/03/2021

Risk of bias in cluster randomised trials

In cluster randomised trials, clusters of subjects such as hospitals or family practices are randomised rather than subjects themselves. [START_REF] Donner | Design and analysis of cluster randomization trials in health research[END_REF] Cluster randomised trials are used for evaluating health service organisation and health policy, often with complex interventions targeted at the level of the cluster, the individual or both. Randomisation should prevent allocation bias at the cluster level provided that it is properly conducted, but differences in individual-level characteristics between the intervention arms can be reintroduced because of the relative timing of participant identification and/or recruitment and cluster randomisation. Indeed, the usual chronology of an individual randomised trial with first recruitment and then randomisation of subjects can be reversed in cluster randomised trials: the identification and/or recruitment of participants frequently take place after randomisation, which could lead to identification/recruitment bias (hereafter called recruitment bias). 2 3 Because blinding is rarely possible for interventions assessed in cluster randomised trials, previous knowledge of the allocation from recruiters and/or participants can influence who is approached and who agrees to participate in a trial. This may lead to different recruitment rates between arms as well as imbalance in participant characteristics. [START_REF] Campbell | Consort 2010 statement: extension to cluster randomised trials[END_REF][START_REF] Hahn | Methodological bias in cluster randomised trials[END_REF][START_REF] Giraudeau | Preventing bias in cluster randomised trials[END_REF] Some solutions proposed to prevent recruitment bias include the identification and recruitment of participants before cluster randomisation or recruitment of participants by a blinded and independent person. [START_REF] Eldridge | Bias in identifying and recruiting participants in cluster randomised trials: what can be done?[END_REF] These solutions should be considered whenever possible, but they are not always feasible or applied. Furthermore, cluster randomised trials are prone to other biases usually encountered when blinding is lacking, namely, performance bias and detection bias. [START_REF] Jüni | Assessing the quality of controlled clinical trials[END_REF] Performance bias refers to systematic differences between the care delivered to experimental and control arms other than the intervention under investigation. Knowledge of the allocation by participants, intervention providers or other trial staff exposes to a risk of performance bias by contamination (with delivery of one of the trial interventions in clusters or participants allocated to receive the other intervention) or by difference in co-intervention delivery. Because they are pragmatic trials, the control arm of cluster randomised trials often consists of usual care (or no intervention). Thus a particular attention may be paid to information provided to participants and care providers in this group in that precise information about the experimental intervention could change their behaviour during the trial and lead to contamination.

Detection bias refers to systematic differences between arms in how outcomes are assessed.

As in individual randomised trials, in cluster randomised trials, knowledge of the allocation by outcome assessors can affect outcome measurement, in particular when outcome 31/03/2021 measurement involves some judgement from the assessor or is directly reported by the participant.

To assess the risk of bias, an accurate description of the distinct procedures is required, but despite existing recommendations [START_REF] Campbell | Consort 2010 statement: extension to cluster randomised trials[END_REF] , the reporting of both the recruitment process and the blinding status for participants and trial staff is often incomplete. 9 10 To help researchers with this description, we developed a graphical tool depicting the sequence and blinding of the different steps of a cluster randomised trial and whether the intervention arms are treated the same or not.

Development of the Timeline cluster tool

The working group consisted of AC, SK, CL, SE, all statisticians, who have been involved in the planning, analysis and reporting of cluster randomised trials as well as in methodological and statistical research on this design. Early in 2015, AC initiated a first version of the graphical tool from real cluster randomised trials. In August of 2015, the working group attended a one-day meeting to discuss this first version of the graphical tool. During the meeting, decisions were made by informal consensus regarding stages that need to be reported and how to better represent cluster and participant levels as well as about blinding. After this first meeting, AC developed a second version of the graphical tool. ET helped AC in the design of this second version (and later versions). Documents were shared by email, and several email iterations took place. Feedback was requested from the whole working group.

We also incorporated feedback received from a presentation of the graphical tool at a meeting on current developments in cluster randomised trials and stepped wedge designs held in London under the auspices of the Royal Statistical Society in October 2015 (all members of the working group were part of the 31 attendees at this meeting) and from pilot use by chief investigators (n=3) and statisticians (n=2) on published or ongoing cluster trials. The latest version of the Timeline cluster tool also takes into account editors' and reviewers' comments.

The Timeline cluster tool

The Timeline cluster tool consists of a diagram and table displayed side by side (see Figures 1, 2 and3 for examples). The diagram represents the sequence of stages of the trial process with successive boxes. Randomisation of clusters is a key stage and is symbolised by a twoway black arrow. All the following stages should be reported when applicable: the identification and recruitment of both clusters and participants, randomisation, intervention 31/03/2021 delivery, and baseline and outcome assessment of participants. The cluster level is symbolised by a ring and the participant level by a stick figure. Blinding status is depicted by the background colour of the boxes. If blinding is complete (all involved protagonists are blinded to allocation), the background of the box is black. If blinding is lacking (no involved protagonist is blinded), the background is white. If blinding is partial (blinding has been used to avoid bias but not all protagonists are blinded to allocation), the background is grey.

Examples of partial blinding include when some protagonists are blinded but others are not or when some protagonists are masked to the hypotheses of the trial. Also stages occurring before randomisation must have a black background not because they are blinded (as no allocation has been made) but because these steps cannot be affected by subsequent allocation. When the control arm receives usual care only, neither symbols appears on the right of the box because nothing is added by the trial at this stage for clusters or participants as compared to standard care. Stages that differ between the intervention arms (at least intervention delivery stage) should be represented by two separate rectangles drawn on each side of the dotted line. The table should at least provide justification for blinding status and other essential details to interpret the diagram, such as the information provided to participants within each arm. The remaining information added in the table is at the user's discretion and replaces what would have been reported in the full text. A Word template to apply the graphical method is provided in the web appendix.

Examples A cluster randomised trial with identification and recruitment of participants before randomisation

The PEACH trial 11 12 is a cluster randomised trial of practice nurses coaching patients with type 2 diabetes. The Timeline cluster diagram for this trial (Figure 1) shows that first clusters are identified and recruited, then participants are identified and recruited, and baseline characteristics are collected before clusters are randomised. All these stages are performed before randomisation, so the corresponding rectangles have black backgrounds and there is no risk of recruitment bias. For the intervention delivery, in the experimental arm, the rectangle has a ring plus a stick figure to the right and in the control arm, neither symbol. We can conclude that the intervention in the experimental arm is delivered at both the cluster and individual levels and consists of usual care in the control arm. The background of the rectangle is white because of no blinding at that stage. The table confirms no blinding for 31/03/2021 general practitioners, practice nurses and participants, so we cannot exclude performance bias.

Finally, blinding is complete for the outcome, change in HbA1c level from baseline to 18 months; there is no risk of detection bias. The equivalent description is about 700 words in the protocol publication [START_REF] Young | Patient Engagement and Coaching for Health: The PEACH study--a cluster randomised controlled trial using the telephone to coach people with type 2 diabetes to engage with their GPs to improve diabetes care: a study protocol[END_REF] and 600 words in the trial report. [START_REF] Blackberry | Effectiveness of general practice based, practice nurse led telephone coaching on glycaemic control of type 2 diabetes: the Patient Engagement and Coaching for Health (PEACH) pragmatic cluster randomised controlled trial[END_REF] A cluster randomised trial with identification and recruitment of participants after randomisation Figure 2 is the Timeline cluster diagram for a cluster randomised trial evaluating a hip protector to reduce hip fractures in older adults. [START_REF] Kannus | Prevention of hip fracture in elderly people with use of a hip protector[END_REF] Before randomisation, clusters are identified and recruited, then participants are identified and assessed. After randomisation, participants are recruited without blinding because recruiters and participants are aware of the allocation and there is a risk of recruitment bias. Blinding is lacking for the intervention delivery targeted at the individual level in the experimental arm and consisting of usual care in the control arm, leading to possible performance bias. Blinding is complete for the primary outcome, hip fracture recorded at the participant level and documented by radiographs. The potential for recruitment bias is confirmed by the trial results, with a significantly higher rate of consent and a lower proportion of participants with a history of falls or severe cognitive impairment at baseline in the control arm as compared to the experimental arm. [START_REF] Kannus | Prevention of hip fracture in elderly people with use of a hip protector[END_REF]

A cluster randomised trial with identification and recruitment of participants after randomisation with measures to prevent bias

The ELECTRA trial 14 is a cluster randomised trial of a specialist nurse intervention to reduce unscheduled asthma care in a multi-ethnic area. Figure 3 displays that even though the identification, recruitment and baseline assessment of individual participants are performed after randomisation of clusters, rectangles for these stages all have a black background because some measures are used to obtain complete blinding and thus prevent recruitment bias. Further details provided in the table indicate that a blinded researcher is used to identify and recruit participants. There is no blinding for the intervention targeted first at the cluster level before participant recruitment, then at the participant level after participant recruitment in both study arms. Some measures are used to avoid detection bias: general practitioners who complete patient records are not blinded, but researchers who extract the primary outcome from the general practice records are blinded, this partial blinding for outcome assessment is represented by the rectangle's grey background.

Comparison with other graphical representations

The Timeline cluster diagram is distinct from the CONSORT flowchart [START_REF] Campbell | Consort 2010 statement: extension to cluster randomised trials[END_REF] which shows the flow of clusters and participants by number approached, randomly assigned, receiving the allocated intervention, and included in the analysis for primary outcome, together with the justification for losses and exclusions. The CONSORT flowchart does not provide information on the chronology of the different stages or blinding of stages. In brief, the CONSORT flowchart is the "How many" and "Why (some participants are excluded)" of the trial, whereas the Timeline cluster diagram is the "When" and "How". Our graphical tool is also distinct from other proposed diagrams that aim to better describe complex interventions either to clarify the timing and differences between arms of their different components (PaT plot method) [START_REF] Perera | Graphical method for depicting randomised trials of complex interventions[END_REF] or to depict interactions between intervention providers at several levels (cascade diagram). [START_REF] Hooper | Cascade diagrams for depicting complex interventions in randomised trials[END_REF] These diagrams provide more detail about the intervention delivery stage only and are useful to enable reproducibility of tested interventions. All these graphical methods are complementary.

Flexibility of the Timeline cluster tool

Our tool offers great flexibility for most cluster trial designs because stages can be added, repeated or removed. For a cluster trial with no need to recruit participants because individual data are obtained from routinely collected data, as in the IRIS trial [START_REF] Feder | Identification and Referral to Improve Safety (IRIS) of women experiencing domestic violence with a primary care training and support programme: a cluster randomised controlled trial[END_REF] , the participant recruitment stage can be removed, with only the participant identification stage retained. Also, a given stage can be performed at different times depending on the intervention arm as in the cluster trial of Cheyne at al., [START_REF] Cheyne | Effects of algorithm for diagnosis of active labour: cluster randomised trial[END_REF] which evaluated the use of an algorithm to diagnose active labour: women consented at admission in the experimental arm and in the postnatal wards in the control arm. This situation of a differential timing for participant recruitment could be represented by a participant recruitment box before intervention delivery in the experimental arm and after intervention delivery in the control arm. For a cluster trial with a repeated cross-sectional design as the one reported by Murphy et al. [START_REF] Murphy | Free healthy breakfasts in primary schools: a cluster randomised controlled trial of a policy intervention in Wales, UK[END_REF]

Discussion

We have proposed a simple and adaptable graphical approach to represent the chronology, blinding and differences between arms in a cluster randomised trial that allows for a quick overview of a given study. The Timeline cluster diagram can be useful at both the design and reporting stages of a cluster randomised trial. At the design stage, the tool may help researchers identify threats to internal validity and consider ways to improve the methodology of their trial, such as use of a recruiter blinded to allocation status. It could also help adequately implementing the trial process in each cluster. At the reporting stage, a more detailed version of the tool can be provided by completing the table with what actually happened during the trial. By providing a precise and adequate description of what was done, it will help readers understand the timeline of trial and appraise the risk of bias. We recommend that future users provide an interpretation of the Timeline cluster diagram in the Discussion section of their trial report. Most often, risk of bias depends on the amount of boxes with black, grey and white backgrounds; indeed, the more black backgrounds, the lower the risk of bias, and the more white backgrounds, the greater the risk of bias, with grey reflecting intermediate or possible risk of bias. However, background colour is not completely associated with the risk of bias level: if the background colour is black, there is no risk of bias at the corresponding stage, but if the background colour is grey or white, the risk of bias must be assessed in the light of other information. For example, if the primary outcome is survival and outcome assessors of vital status are not blinded to allocation, the box for outcome assessment must have a white background even though there is only low risk of detection bias.

Weaknesses in design and conduct as well as incomplete reporting of biomedical studies are important leading causes of research waste. 21 22 Risk of bias, specifically recruitment bias, is an issue in cluster randomised trials, with remaining room for improvement. We believe that our graphical approach could help achieve a better management and reporting of cluster randomised trials, allowing for an informed assessment of the risk of bias. We have received positive feedback from the investigators and statisticians who have used the current version of the Timeline cluster tool, but we anticipate that some further enhancement will probably be 31/03/2021 required. Therefore, we encourage suggestions from readers and feedback from the practical experience of future users.

  in which clusters are followed over time but participants change during the trial, the need to repeat the identification and recruitment of participants after intervention delivery could be represented by a loop starting and ending on the central dotted line of the diagram. A cluster cross-over trial such as REPHVIM trial 20 could be depicted by adding two crossing and ascending arrows to depict the switch from one intervention to the other at the end of the first period. 31/03/2021 Timeline cluster diagrams corresponding to these aforementioned examples are provided in the web appendix.

Acknowledgements We thank Catherine Laporte, Chris Gale and Ken Haguenoer, investigators of cluster randomised trials, Patty Chondros and Karla Hemming, statisticians, for testing and giving feedback on the Timeline cluster graphical tool. We also thank the reviewers of this manuscript for helpful comments.

Contributors and sources: BG had the idea of the graphical tool. AC conceived the graphical tool and drafted the original manuscript. SK, ET, CL, SE and BG helped in the creation of the graphical tool and contributed to the final version of the manuscript. AC is the guarantor.

Competing interests:

We have read and understood the BMJ Group policy on declaration of interests and declare we have no competing interests. Data sharing: A Word template is available through the BMJ website or on request from the corresponding author I Agnès Caille, the Corresponding Author has the right to grant on behalf of all authors and does grant on behalf of all authors, an exclusive licence on a worldwide basis to the BMJ Publishing Group Ltd to permit this article (if accepted) to be published in BMJ editions and any other BMJPGL products and sublicences such use and exploit all subsidiary rights, as set out in our licence.