François Laurent
email: francois.laurent@pasteur.fr

Charlotte Floderer

Cyril Favard

Delphine Muriaux

Jean-Baptiste Masson
email: jean-baptiste.masson@pasteur.fr

Christian L Vestergaard
email: christian.vestergaard@pasteur.fr

Mapping spatio-temporal dynamics of single biomolecules in living cells

est destinée au dépôt

INTRODUCTION

Photoactivation fluorescence microscopy makes it possible to track and study the dynamics of single biomolecules in living cells and at the whole cell scale [1,2]. Such single molecule tracking approaches have revealed rich intracellular transport dynamics [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. In many cellular systems the molecules' dynamics show spatial heterogeneities which have been linked to biological function. Recent examples include neuroreceptors binding in the postsynaptic membrane [9-11, 13, 16], signaling molecules involved in the initiation of lamellopodium formation in motile cells [18], protein transport in the endoplasmic reticulum [19], and recruitment of polyproteins in the plasma membrane to form HIV-1 virus-like particles [20]. A successful approach to analyse single molecule tracking data from such systems is to assume that each tracer, of the same molecular species, experiences the same spatially dependent force/drift and diffusivity fields. This hypothesis makes it possible to infer spatial maps of the tracers' dynamics from many, typically short, recorded trajectories [7,[22][23][24][25][START_REF] Frishman | Learning force fields from stochastic trajectories[END_REF]. Recent experimental modalities can produce datasets consisting of up to tens of millions single molecule localisations. This opens up the possibility to map not only the spatial variation of biomolecular dynamics but also study how these maps evolve in time [20].

When inferring such maps, a fundamental trade-off between spatial (and temporal) resolution and statistical precision results from the limited amount data. On the one hand, the spatial and temporal resolutions should be as high as possible in order to capture the spatial and temporal phenomena of interest. On the other hand, the effective number of parameters should be kept as low as possible in order to be able to reliably infer them from the available data.

The spatio-temporal scale of mapping is typically set by binning the recorded localisations. To apply the mapping approach to processes of a priori unknown spatial and temporal scales and mitigate the impact of arbitrarily choosing such a scale, we propose to regularise the inferred maps. This can be done in a Bayesian setting using field regularising priors [13,24]. It increases the precision and robustness of inferred parameter maps by imposing that they must vary smoothly, and makes it possible to infer high resolution maps even from sparse data.

We here extend the Bayesian mapping framework of [7,13,24] to include both spatial and temporal regularisation for robust and statistically principled inference of time-varying large-scale spatial maps of the diffusivity and potential fields governing single biomolecule motion. The price to pay for the increased precision offered by regularisation is a drastic increase in the computational complexity as it couples the parameters of the map. This transforms a collection of independent low-dimensional inferences of individual map parameters into a global high-dimensional inference problem. The classic approach for parameter inference such as direct gradient based algorithms or Markov chain Monte Carlo (MCMC) are ill suited for high dimensional inferences from large datasets since they scale poorly with the size of the dataset and the number of parameters to infer.

To tackle the computational complexity of inferring large spatio-temporal maps, we instead make use of stochastic optimization [START_REF] Robbins | [END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF], a general yet powerful method for addressing high-dimensional optimization problems on large datasets. Stochastic optimisation can be traced back to 1951 [START_REF] Robbins | [END_REF], but has more recently seen widespread application to handle the massive datasets found in modern statistical and machine learning problems [START_REF] Kingma | Auto-encoding variational bayes[END_REF].

In this article, we describe a stochastic optimisation procedure tailored to the mapping problem, which uses spatio-temporally local mini-batches to increase the efficiency of updates and enable multiprocessing. We validate the approach on numerical data and demonstrate the advantage of stochastic optimisation over direct optimisation of the global posterior.

Although we focus on maximum a posteriori (MAP) inference of map parameters, we also provide a procedure for sampling the posterior distribution based on stochastic gradient Langevin dynamics [START_REF] Welling | Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML'11[END_REF] in the Supplementary Information.

We apply the framework to infer spatio-temporal maps of the intracellular dynamics of Gag proteins, which are involved in the formation the HIV-1 capsids at the plasma membrane of T cells. In this article, we show transient aggregation events that aborted before they led to fully-formed virus-like particles (VLP).

The Bayesian framework for spatio-temporally regularised mapping of single-molecule dynamics developed here is implemented in the open-source and freely available TRamWAy project (github.com/DecBayComp/ TRamWAy).

METHODS

Model

We model the single molecule dynamics using the heterogeneous overdamped Langevin equation (OLE):

dr(t) dt = a(r, t) + 󰁳 2D(r, t) ξ(t) . (1)
Here a(r, t) and D(r, t) are the (spatially and temporally varying) deterministic drift and diffusivity of the molecules in the point r at time t, respectively, and ξ(t) is a continuous-time white noise process, defined as the derivative of the Wiener process.

In the applications presented here, we will assume that the dynamics are effectively potential-driven and that the fluctuation-dissipation relation is respected locally. The first assumption imposes that the drift should be proportional to the gradient of a scalar potential field, a(r, t) = -∇V (r, t)/γ(r, t). The second assumption links the proportionality constant, the drag coefficient γ(r, t), to the diffusivity through the Einstein relation γ(r, t) = k B T /D(r, t).

These assumptions are typically employed to interpret the inferred drift fields as stemming from effective potential energy landscapes [7,9,13,16,20,[START_REF] Frishman | Learning force fields from stochastic trajectories[END_REF][START_REF] Chang | [END_REF]32]. Note however that since equilibrium conditions are not guaranteed in biological systems, inferred potential maps cannot generally be interpreted as potential energies in a strict sense as they may include terms of unknown amplitude due to spatial variations in the diffusivity [START_REF] Serov | Robust inference of forces in heterogeneous environments[END_REF] and terms induced by non-equilibrium energy fluxes. Notwithstanding, the inferred maps may capture biologically relevant information regardless of whether the equilibrium assumption is satisfied or not [9,13,16,20]. As a statistical description, a model that only assumes potential-driven dynamics, and not equilibrium, is equivalent to the equilibrium model. Note also that assumptions about equilibrium can be relaxed significantly by using the approach developed in [START_REF] Serov | Robust inference of forces in heterogeneous environments[END_REF] for non-regularised maps.

Bayesian inference and regularisation

Bayesian inference deals with estimating the posterior probability p(θ|x) of model parameters θ conditioned on data x. Here, the data is a set of M single molecule trajectories,

x ≡ r = {{r m n } Nm n=1 } M m=1
, where r m n denotes the nth recorded position in the mth trajectory, and N m the trajectory's length.

The posterior is related to the model likelihood p(x|θ) through Bayes' theorem,

p(θ|x) = p(x|θ)π(θ) p(x) , (2)
where π(θ) is the prior probability for the parameters θ in absence of data and p(x) = 󰁕 p(x, θ)dθ is a normalising constant, termed the evidence. The evidence p(x) can be ignored for parameter inference but is central to Bayesian model selection [START_REF] Kass | [END_REF].

Equation (2) provides a probabilistic framework for regularisation by employing smoothing priors instead of conjugate or noninformative priors (see below).

Maximising p(θ|x) w.r.t. θ leads to maximum a posteriori (MAP) estimation. In this article, we will focus on MAP estimation, but we also provide a procedure for sampling the posterior around the MAP in the Supplementary Information.

Spatio-temporal meshing and approximate likelihood

The likelihood for our model can be obtained as the fundamental solution (the Green's function) of the Fokker-Planck equation corresponding to Eq. ( 1). However, in general such a solution is analytically intractable and numerically too time-consuming to be considered for inference.

Instead, we follow and extend the approach developed in [7,24] and approximate the fields V (r, t) and D(r, t) as piecewise constant in both t and r. In practice, we tessellate/segment space and time into spatio-temporal domains and consider instead the sets of average values of V (r, t) and D(r, t) in each domain, V = {V α,τ } and D = {D α,τ }, where τ = 1, 2, . . . , T indexes the time segments, and α = 1, 2, . . . , Ω identifies the space domains.

⌧ + 1 ↵ ⌧ ⌧ 1
FIG. 1. Illustration of the spatio-temporal tessellation. We used regular meshes consisting of equal-sized hexagons in space and equal-length time segments. The highlighted domains (mesh domains with thick black edges) show the spatiotemporal neighbourhood Bα,τ of the central domain (α, τ ), with (α, τ ) itself included. We refer to the spatial neighbourhood of a spatial domain α (disregarding time) as Rα, and to the temporal neighbourhood of a temporal segment τ as Tτ = {τ -1, τ, τ + 1}. Temporal neighbours share the same spatial position and spatial neighbours are considered at fixed time.

Figure 1 illustrates the spatio-temporal tessellation we employ here, which consists in a regular hexagonal tessellation combined with a sliding time window. It also illustrates the important concepts of the spatial neighbourhood R α of a space domain α, containing α itself and the space domains in contact with α, the temporal neighbourhood T τ = {τ -1, τ, τ + 1} of a time segment τ , and the spatio-temporal neighbourhood B α,τ = {(α ′ , τ )} α ′ ∈Rα ∪{(α, τ ′ )} τ ′ ∈Tτ of spatio-temporal domain (α, τ ). Note that the domains can overlap in space and time (see Data sampling below).

Many other geometries are possible for the domains, which need not be regular and can be adapted to the density of localisations. We here consider only regular hexagonal tessellations, which have several attractive properties for spatio-temporal mapping. Using a mesh that is constant over time makes it straightforward to implement temporal regularisation as one does not have to deal with partially overlapping domains (Fig. 1). Similarly, the hexagonal spatial meshing avoids having to arbitrarily choose how to make regularisation of a domain depend on the distances to the centers of neighbouring domains and on their shared perimeters since these are all equal. In the case of a square meshing this problem is especially egregious as neighbours along the diagonal meet in a single point only. Instead of using an adaptive tessellation, we rely on regularisation to handle the domains with little data. Indeed, regularising makes the inference more robust to noisy data and allows to consider higher resolutions for the entire maps without having to sacrifice either statistical precision in low density regions or spatio-temporal resolution in high density regions.

Under the piecewise constant approximation, the fundamental solution to the Fokker-Planck equation is a Gaussian distribution, and the probability for a tracer located in the spatio-temporal domain (α, τ ) to perform a displacement of ∆r i is then [7,24]:

p(∆r i |θ Rα,τ ) = 1 4πD α,τ ∆t i exp 󰀳 󰁃 - 󰀏 󰀏 󰀏|∆ri + Dα,τ ∆ti kBT ∇V α,τ 󰀏 󰀏 󰀏 | 2 4D α,τ ∆t i 󰀴 󰁄 .
(3) Here, the potential gradient ∇V α,τ is evaluated numerically, in practice using a finite differences scheme. It thus depends not only on the local parameters θ α,τ = (D α,τ , V α,τ ) in the domain, but on the values of V in the whole spatial neighbourhood of (α, τ ):

θ Rα,τ = {θ α ′ ,τ } α ′ ∈Rα .
The total likelihood of all recorded trajectories can be written in a factorised form as a product of local likelihoods [Eq. ( 3)]. If we let ∆r α,τ = {∆r α,τ,i } i denote the set of all displacements recorded inside the domain {α, τ }, we can write the total likelihood for all the recorded data as:

p(∆r|θ) = 󰁜 α,τ p(∆r α,τ |θ Rα,τ ) , (4)
where θ = {θ α,τ } is the set of all local map parameters and ∆r = {∆r α,τ } is the set of all recorded displacements, with α = 1, 2, . . . , Ω and τ = 1, 2, . . . , T . Note that for the factorised likelihood [Eq. ( 4)], ∆r is equivalent to the set of trajectories r = {{r m n } Nm n=1 } M m=1 .

Spatio-temporal regularisation

To regularise the inferred maps, we impose a field regularising (Sobolev) prior π(D, V ) that penalizes variations in D and V between neighbouring domains in space and time [35,[START_REF] Engl | Regularization of inverse problems[END_REF]. This prevents physically-unrealistic patterns of field parameters such as spurious minima or maxima. The regularisation prior factorizes as π(D, V ) = π(D)π(V ), with

π(D) ∝ e - 󰁕 (µr||∇D|| 2 +µt Ḋ2 )dt dr , (5)
π(V ) ∝ e - 󰁕 (λr||∇V || 2 +λt V 2 )dt dr . (6)
In practice, since the maps are modeled as being piecewise constant, the integrals above reduce to sums:

π(D) ∝ e - 󰁓 α,τ δτ Aα(µr||∇Dα,τ || 2 +µt Ḋ2 α,τ ) , (7)
π(V) ∝ e - 󰁓 α,τ δτ Aα(µr||∇Vα,τ || 2 +µt V 2 α,τ ) . ( 8
)
where δ τ is the temporal duration and A α is the area (in 2D) or volume (in 3D) of domain (α, τ ). For regular maps, δ τ = δ and A α = A are constant, but we allow them to vary in general for applications with non-regular meshing.

In π(D) and π(V), the differential operators ∇ and involve calculating finite differences [see Stochastic optimisation, Eqs. ( 12) and ( 13), below].

Including the regularising priors, the posterior reads:

p(θ|∆r) ∝ p(∆r|θ)π(D)π(V) . (9)
As seen from Eq. ( 4), (7), and ( 8) neither the likelihood nor the priors factorise into independent local likelihoods. This couples the inferences of map parameters of the different domains, both spatially and temporally, which in turn renders the optimisation of the posterior computationally hard. Indeed, we cannot simply optimise the posterior p(θ|∆r) w.r.t. the parameters in each map domain independently, but instead need to optimise w.r.t. all parameters simultaneously.

Stochastic optimisation

Here, we propose to rely on stochastic optimisation [START_REF] Robbins | [END_REF][START_REF] Spall | Introduction to Stochastic Search and Optimization[END_REF] to solve our inference problem. In stochastic optimisation, the cost function-here the negative log posterior f (θ) = -log p(θ|∆r)-is optimised by following a set of noisy estimates of its gradient based on minibatches of data.

In its general form, stochastic optimisation randomly selects a subset of n ≪ N points, out of the total dataset of N points, to calculate the gradient of the cost function at each optimisation step. This can dramatically decrease the computational burden of high dimensional inferences from massive datasets.

We here describe our algorithm briefly and refer to the Supplementary Information for technical details and implementation procedures. We devised a stochastic optimisation scheme to take advantage of the fact that the parameters of each domain only couple to spatio-temporally neighbouring domains. This lets us decompose the cost function f as a sum of local terms,

f (θ) = 󰁛 α,τ f α,τ (θ Bα,τ ) , (10)
where each local term depends on the neighbourhood

B α,τ = {(α ′ , τ )} α ′ ∈Rα ∪ {(α, τ ′ )} τ ′ ∈Tτ , comprising (α, τ )
as well as spatio-temporally neighbouring domains (see Figure 1). Each local cost term is given by:

f α,τ (θ Bα,τ ) = -log p(∆r α,τ |θ Rα,τ ) (11) +µ r q α (D Rα,τ ) + µ t q τ (D α,Tτ ) +λ r q α (V Rα,τ ) + λ t q τ (V α,Tτ ) .
Here, the local regularisation terms q α and q τ are given by:

q α (φ Rα,τ ) = 1 2 δ τ A α |R α | -1 󰁛 α ′ ∈Rα󰄀{α} (φ α,τ -φ α ′ ,τ ) 2 ||r α -r ′ α || 2 , ( 12
)
where φ is a placeholder for either D or V and r α and r ′ α are the centers of the space domains α and α ′ , respectively, and

q τ (φ α,Tτ ) = 1 2 δ τ A α |T τ | -1 󰁛 τ ′ ∈Tτ 󰄀{τ } (φ α,τ -φ α,τ ′ ) 2 ∆T 2 . ( 13
)
The factor 1/2 in Eqs. ( 12) and ( 13) accounts for the fact that each individual penalty term is counted twice [e.g.

(φ α,τ -φ α,τ ′ ) 2 /∆T 2 is counted in both f α,τ and f α,τ ′ ].
The local cost functions f α,τ are not independent, but the decomposition [Eq. (10)] enables partial evaluation of f for gradient calculation and line search in parameter subspaces. We devised a stochastic algorithm, based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization method [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF], that only requires such partial evaluations of f . This algorithm updates a subset of the parameters based on a local minibatch of data at each iteration. Rather than randomly sampling data minibatches, the algorithm selects data that are directly related to a set of local physical parameters. This selection is thus only local and directly maps onto the mesh domains associated with the parameters. The locality of the minibatches also makes the algorithm amenable to parallelisation and multiprocessing (see Supplementary Information for details on multi-process implementation).

The optimisation algorithm repeats the following steps until convergence: For each iteration k, 1. select a local domain (α, τ ) to update; 2. sample a minibatch ∆r Bα,τ = ∪ (α ′ ,τ ′ )∈Bα,τ ∆r α ′ ,τ ′ of displacements in (α, τ )'s neighbourhood;

3. update the local parameters following θ

(k) α,τ = θ (k-1) α,τ + ∆θ(∆r Bα,τ , θ (k-1)
Sα,τ ), where ∆θ is an implicit function of the local data ∆r Bα,τ and the parameters in the extended neighbourhood

S α,τ = ∪ α ′ ,τ ′ ∈Bα,τ B α ′ ,τ ′ (see Gradient descent in a sub- space below).
The order in which the domains are updated in (1) should be randomised in order to avoid artefacts in the inferred maps. We do so by looping over all domains in an order that is randomised in each epoch. This ensures that all domains are updated while avoiding an arbitrary ordering of the updates. We here select the whole set ∆r Bα,τ as minibatch in (2), but it is possible to subsample it instead for large datasets. The update ∆θ of θ α,τ in (3) only depends on the local cost terms f α,τ in the neighbourhood B α,τ . Since each local cost function f α ′ ,τ ′ depends on parameters of the domains in its neighbourhood B α ′ ,τ ′ , the update ∆θ ultimately depends on the parameters of the extended neighbourhood S α,τ that includes all neighbours of the domains in B α,τ , but not any domains further away.

Gradient descent in a subspace

A detailed description of how ∆θ is calculated is given in the Supplementary Information. Briefly, for each iteration k we obtain ∆θ from the implicit relation ∆θ(∆r Bα,τ , θ

(k-1) Sα,τ ) = s (k) p (k) . Here s (k) is the step size and p (k) = -H (k-1) ∇ α,τ f (θ (k-1)
) is the descent direction, with ∇ α,τ f the gradient of f w.r.t the local parameters θ α,τ and H an approximation to the inverse Hessian matrix. We relied on approximate line search [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] to estimate s (k) . Interestingly, the same partial evaluation as for gradient calculation stands, which makes line search operate fast, in a 2-dimensional parameter subspace.

Data sampling

This study showcases three datasets of simulated trajectories and one experimental dataset consisting of trajectories of wild-type Gag proteins (WT Gag) recorded in the plasma membrane using sptPALM. In this section, we briefly describe how these datasets were binned and how molecule locations were assigned to the resulting spatiotemporal domains.

For each dataset, a regular grid of non-overlapping hexagons was adjusted to the bounding box of all the molecule locations. The minimal diameter (the diameter of the inscribed circle) of the hexagons was set to 100 nm for the first two simulations, 68.5 nm for the third simulation, and 20 nm for the WT Gag dataset. This diameter defines the spatial resolution of the analysis.

Similarly, the total recording (or simulation) time was divided into regular intervals. The first simulated dataset did not feature time-evolving dynamics, so only a single time-slice was used. The second and third simulated datasets feature 18 time segments of 1 s and 5 s respectively, while the experimental WT Gag dataset features 12 1-min time segments. The experimental modalities of the WT Gag dataset and statistical analysis is described in more detail in the following section (Experimental data).

To further ensure the reliability of the inference procedure for sparse data, we incorporated different oversampling modalities that consider overlapping domains and may consequently assign localisations to multiple domains. Namely, (i) considering any initial tessellation, the spatial domains were enlarged as overlapping circular regions of constant radius; (ii) the temporal extent of the domains was defined using a sliding time window that introduced some constant overlap between successive time segments; and (iii) the domains that did not include a given minimum number of localisations were individually and symmetrically enlarged in time until this minimum number was reached. While the oversampling has a smoothing effect, the main purpose is to ensure the reliability of the optimisation and inference procedures. For visualisation, the non-overlapping shapes from the original tessellation were shown. In all three oversampling modalities, the spatial and temporal integration in q α and q τ (respectively) considered the non-overlapping extent of the domains, e.g. for q α in modality (i) the original non-overlapping shapes of the spatial domains.

No oversampling was performed for the first simulated dataset. For the second simulated dataset, following modality (iii) spatio-temporal domains with less than 20 molecule locations were individually enlarged in time until 20 locations could be assigned to the domain. No oversampling was performed for the third simulated dataset in order to isolate the effect of the smoothing priors. For the experimental dataset, regular spatial (i) and temporal (ii) windowing were applied. Since the dataset is very sparse and our focus is to detect short-lived events, we used relatively large spatial overlap (130 nm for 20 nm spatial resolution) and little temporal overlap (2 min sliding windows for a 1 min temporal resolution), and we did not apply adaptive temporal windowing.

Note that while the regularizing priors in principle suffice to handle domains with few datapoints, it is in practice preferable to exclude domains not containing a minimal number of localisations and/or artificially increasing the number of points by oversampling in order to ensure that inferences are stable and the posterior is integrable. This is because each spatio-temporal domain contributes equally to the regularisation term in the posterior probability no matter how many points it contains. Thus, if many domains contain few or no points (as is the case in the experimental dataset), these will contribute disproportionally to the posterior, and in a manner that is overly sensitive to initialisation, possibly rendering the inference procedure unstable.

Experimental data

We analyse a previously unpublished set of trajectories of wild-type (WT) mEOS2-tagged HIV-1 Gag proteins recorded at 50 Hz in the plasma membrane of a single CD4 + T cell for a previous study [20]. The experimental procedure, from protein production (using the pGag-(i)mEOS2 WT plasmid) to imaging and tracking, is described in [20] (condition: WT Gag).

To map the single Gag dynamics, here the spatial 2D domains were regularly laid across a 3 µm × 3 µm area, with centers arranged along a hexagonal grid, spaced by 20 nm and sized as 130 nm-radius discs. This large radius-as large as observed virus-like particles (VLPs) [20]-was necessary to collect enough molecule locations so that the diffusivity and potential energy could be estimated over dense and large spatio-temporal patches that cover not only patterns of interest, but also the close environment of these patterns. To emphasize the benefit of temporal regularisation, we used a 2 min sliding temporal windows and limited the overlap between successive time segments to 50%, resulting in one time segment every minute. Spatio-temporal domains containing no localisations were excluded from the analysis.

The effects of localization uncertainty and motion blur (also known as static and dynamic errors, respectively [START_REF] Savin | [END_REF]) were accounted for by replacing D(r, t) by D eff (r, t) = (1 -2R)D(r, t) + σ 2 /∆t in Eq. ( 3) [39]. Here σ 2 is the variance of the localization uncertainty, assumed zero mean and time-uncorrelated, and R is a coefficient quantifying the motion blur [40,41], where R = 1/6 corresponds to continuously open camera shutter and R = 0 may be obtained by leaving the camera shutter open for only a vanishing part of the time-lapse between measurements. For the present measurements, the localisation uncertainty was estimated to be σ = 30 nm and the motion blur coefficient was R = 1/6.

RESULTS

Stochastic optimisation accelerates inference of spatio-temporal maps

We first demonstrate the stochastic optimisation approach on simulated single particle tracking data, and we compare it to the deterministic global optimisation. We namely consider a low-dimensional example of inference of a static, spatially varying potential field (558 parameters to infer and 2 682 particle localisations) and a higherdimensional intermediary example of a spatio-temporally varying field (4 648 parameters, 22 589 particle localisations). Both datasets exhibit constant diffusivity and a potential energy sink in space, which is constant in time for the first dataset (Fig. 2) and varies over time for the second (Fig. 3). The first dataset was introduced to quantify the convergence of several optimisation methods and does not involve time windowing. The second dataset introduces time regularisation and emphasizes the additional benefit of stochastic optimization in this context. It is divided into 18 non-overlapping time segments as illustrated in Fig. 3.

Convergence is measured by the correlation of the partial estimates φ (k) (where φ ∈ {D, V}) of each method with the final estimate φ * obtained using the global deterministic optimisation of the full cost function. This way we compare not only the speed of convergence of the methods, but also that they converge towards the same optimum.

Fitting V took more local cost function calls than fitting D did (Figs. 4 and5). This is generally the case even when D exhibits equivalent patterns. We will therefore focus on convergence for V.

Spatial map

The lower computational complexity of all the improved methods compared to the naive gradient calculation is illustrated in Figure 4. Direct gradient descent converges much slower than its quasi-Newton counterpart, although this variant can also run in multi-processing mode and consequently be faster (in absolute processing time) than the non-stochastic variant. This validates the approximation made in estimating the inverse Hessian matrix as a diagonal 2x2-block matrix (see Supplementary Information).

All optimisation methods converge relatively fast for the inference of the low-dimensional spatial map. The stochastic variant initially converges fastest until some precision is reached (for the potential energy V at a correlation of 0.9997 with the optimal parameter values). Above this level, the non-stochastic variant converges at a constant pace whereas the stochastic variant slows down, perhaps due to many local updates being redundant when most parameters are close to their optimal values. This is mitigated by the parallelisability of stochastic optimisation. Since multiprocessing does not increase the number of local cost function calls (Figure 4), this allows to divide the elapsed time for optimisation by the number of workers (about 20 on modern computers) using multi processing as compared to single processing.

Spatio-temporal map

The second dataset involves a transient potential energy sink as illustrated in Figure 3. The spatial-temporal map we aim to infer contains 18 individual time slices, so the number of parameters is 18 times higher than for the static map in the low-dimensional example above.

We compared stochastic and non-stochastic quasi-Newton optimisation (Figure 5). For the regularised spatio-temporal inference, the stochastic algorithm converges much faster than the non-stochastic one-about an order of magnitude for the present example-and the parallelisability of the stochastic algorithm compounds this gap by another order of magnitude.

Note that the number of parameters inferred in this example (4 648) is much smaller than the number of FIG. 3. Inferred potential energy landscapes at successive time segments in the intermediary spatio-temporally varying example (4 648 parameters in total). FIG. 4. Convergence as function of the number of calls to local cost functions fα,τ for the low-dimensional spatial inference (558 parameters to infer). n.s. naive stands for non-stochastic with naive gradient calculation. H = I means that the inverse Hessian matrix is not approximated using the BFGS update rule, but set to the identity matrix instead; this is the direct gradient descent, in contrast to the quasi-Newton approach. s.p. and m.p. stand for single-process and multi-process, respectively. FIG. 5. Convergence as a function of the number of calls to local cost functions fα,τ for the spatio-temporal inference (4 648 parameters to infer). m.p. stand for multi-process.

parameters one may meet in empirical examples. The map inferred in the experimental example presented below for example contains 127 578 parameters. For such high-dimensional inferences, the non-stochastic algorithm does not converge.

Computational complexity

We let P denote the number of parameters. P is proportional to the number of different domains (α, τ ) in the mesh, so the full gradient calculation of non-stochastic algorithms has a computational complexity of O(P 2 ) per iteration. In comparison the partial gradient evaluation of our stochastic algorithm has a complexity of O(P B), where B = |B α,τ | is the size of a local neighbourhood. In terms of the size of the parameter space, P , this reduces to O(P ) since B is always bounded by a low value in practice: here B = 9 and typical values in 2D are B < 10. For reference, P = 127 578 in the experimental example below.

Strategies for choosing regularisation hyperparameters

Figure 6 illustrates how the value of the hyperparameter λ τ controlling the strength of temporal regularisation of the potential field V influences the inferred field. The dataset involved is slightly larger then the previous one. It features 9 992 parameters and 126 022 localisations, and it displays three potential wells (Fig. 6, top rows).

The optimal choice of regularisation hyperparameters µ r , µ t , λ r and λ t depends on the overall study, and there is no general be-all and end-ell strategy for choosing them. The application of mapping to biomolecule dynamics in living cells may involve the characterisation of specific features of the inferred maps, such as the size, depth and duration of potential wells. Some studies then focus on differences between various biological conditions rather than on absolute measurements. More generally, inferring maps is often a preliminary analysis step, and the further analysis steps may all bring their specific requirements, which can guide the selection of regularising hyperparameters. A common criterion is to avoid false detections of differences between the feature values corresponding to different biological conditions. False positives may notably be caused by unphysical outliers in the distributions of values of map features as statistical tests used for comparison often are sensitive to outliers. A simple approach to correct this consists in tuning the hyperparameters to ensure that the resulting distributions do not exhibit long tails, as measured for example using a Kolomogorov-Smirnov test to compare the empirical distribution with a Gaussian fit [20]. As an undesired consequence of using smoothing priors, the inferred parameter maps may exhibit biased (flattened) values.

Conversely, if we want to detect sharp patterns or extreme values, the hyperparameters should be tuned to the smallest values that still let the inference converge and make the inferred maps exhibit as few holes (undefined parameter values) as possible within the regions of interest.

Another approach, applied in [13] (Supplementary Information), tunes the spatial regularisation hyperparameters using large-scale numerical simulations based on realistic parameter landscapes. Numerical trajectories are generated so that they mimick experimental trajectories, and a grid search is performed to find the hyperparameter values that minimise the discrepancy (reconstruction error) between the inferred and expected maps.

Multiple procedures for tuning hyper-parameters can be envisioned depending on the nature of the experiments and the expected results. While we recommend simulation-based exploration of the effects of hyperparameters, increasing the robustness of statistical tests is often taken as a priority in the analysis of biological data. Generally the two procedures are not mutually exclusive, but simulation-based tuning can be computerintensive. Furthermore, if the regions of biological interest are highly local in space and/or time, minimising the reconstruction error will tend to put most weight on domains that are of lesser interest because there are more of these.

The experimental dataset introduced in this article is an example of such a complex inference with sparse regions of interest. The hyperparameters were set to the following values: µ r = 1, µ t = 10, λ r = 0.1 and λ t = 1.

Spatio-temporal mapping of Gag protein dynamics

We illustrate our inference framework by inferring a spatio-temporal map of the dynamics wild-type mEOS2tagged HIV-1 Gag proteins in a single CD4 + T cell. The Gag protein is involved in the budding of HIV-1 virus particles at the plasma membrane and, when expressed alone, is sufficient to form inactive virus-like particles (VLPs) [20].

We mapped the dynamics in a 3 µm × 3 µm square of the plasma membrane over 15 min at spatial and temporal resolutions of 20 nm and 1 min (see Methods for a detailed description of the dataset and analysis). To observe transient interactions with possibly shorter du- rations than the classically observed 20-minute particle budding, we segmented time using a 2-minute sliding time window and a 1-minute resolution. This is to be compared to previous analyses that ran a 4-minute time window [20].

The dataset exhibited low density-white areas in Fig. 7 were not visited by any molecule during the corresponding time segment. Clear patterns could still be identified, such as the ∼ 100 nm energy trap illustrated on the top row. This trap lasted less than 2 minutes, which is indicative of transient assembly that did not lead to a fully assembled VLP [20]. It is worth noting that these short events cannot be detected based on the local density of Gag and that information still lacks on VLP that did not assemble. Other weaker traps were also identified, such the ones illustrated on the third and fourth rows.

A key benefit of regularising in time is the homogeneity of inferred potential energies between distant locations and distinct time segments. This is not the case when individually inferring spatial maps for different time segments since potential energies are only determined up to an additive constant. The temporal regularisation made these estimated energies directly comparable between one another.

In spite of the 127 578 parameters, the inference took less than 24 hours on a desktop computer with an Intel c CONCLUSIONS Here, we introduced a Bayesian framework for fully regularised inference of spatio-temporal maps of intracellular biomolecule dynamics from their recorded trajectories. To solve the computational optimisation problem of inferring the most probable (maximum a posteriori -MAP) values of the thousands of map parameters we relied on stochastic optimisation -a general class of optimisation methods well suited to high-dimensional inferences and large datasets. We exploited the particular physical structure of the mapping problem by adapting the quasi-Newton BFGS algorithm to make use of spatio-temporally local minibatches, which provide particularly efficient updates as they involve only small parameter subspaces while separating experimental localisations that are linked to a limited set of parameters. While our focus was on MAP estimation, we additionally provided a complementary procedure based on a modern technique from Bayesian statistics and machine learning to access the full posterior distribution: sampling of the posterior around the MAP using stochastic Gradient Langevin dynamics. We demonstrated the advantage of stochastic optimisation over direct optimisation of the full posterior, and, for the first time, we mapped the spatio-temporal dynamics of HIV-1 Gag proteins on large areas with temporal regularisation.

θ (k) = θ (k-1) + 󰂃 2 ∇P (θ|{∆r}) + √ 󰂃 η k , (14)
where η k follows a standard Gaussian distribution, η k ∼ N (0, 1). Equation ( 14) can be recognised as the Euler approximation to an overdamped Langevin equation with the posterior distribution acting as an effective potential and 󰂃 as a diffusivity. This procedure was shown to allow efficient sampling of the posterior distribution, and it converges quickly when initialised at the MAP values of the parameters. Note that since we here initialise the procedure at the MAP values for θ, we do not need to make 󰂃 vary with time as in the original publication [START_REF] Welling | Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML'11[END_REF] where the scheme was also used to optimise the posterior before sampling around the MAP.

STOCHASTIC OPTIMISATION

We here give a detailed description of the proposed stochastic optimisation algorithm. Each physical parameter φ (indexed by ℓ) in the set of model parameters θ, namely the diffusivity D and the potential energy V , is represented by as many model parameters as spatial subdomains (indexed by α) times temporal segments (indexed by τ ). We may write the set of parameters formally as:

θ = ∪ ℓ θ ℓ = ∪ α,τ ∪ ℓ {θ ℓ,α,τ } = ∪ α,τ {D α,τ , V α,τ }.
To show the generality of the stochastic approach, we here consider a more general formulation that in the main text by letting B α,τ be any neighbourhood of (α, τ ), instead of just being comprised of nearest neighbours, and by letting θ α,τ be any collection of local parameters, e.g. mobility and drift θ α,τ = (ψ α,τ , a α,τ ) or drag and force θ α,τ = (γ α,τ , f α,τ ). In this more general case, the local cost can be written as

f α,τ (θ Bα,τ ) = -log(P (θ Rα,τ |{∆r} α,τ )) (15)
+ 󰁛 ℓ [µ ℓ,r q α (θ ℓ,Rα,τ ) + µ ℓ,t q τ (θ ℓ,α,Tτ )] .

Parsimonious gradient calculation

The cost function f being a sum of local components makes partial evaluations of the gradient possible. For example, considering any increment 󰂃 > 0 along a single scalar parameter θ + = θ + 󰂃 1 ℓ,α,τ and θ -= θ -󰂃 1 ℓ,α,τ , we have:

f (θ + ) -f (θ -) = 󰁛 α ′ ∈Rα 󰁫 f α ′ ,τ 󰀓 θ + B α ′ ,τ 󰀔 -f α ′ ,τ 󰀓 θ - B α ′ ,τ 󰀔󰁬 + 󰁛 τ ′ ∈Tτ 󰁫 f α,τ ′ 󰀓 θ + B α,τ ′ 󰀔 -f α,τ ′ 󰀓 θ - B α,τ ′ 󰀔󰁬 . (16)
Only a few local functions {f α,τ } are evaluated above. The above calculus is the basis for gradient calculation. Indeed, the component (ℓ, α, τ ) of the gradient ∇f is calculated as: ∇f ℓ,α,τ (θ) = f (θ + )-f (θ -) 󰂃 1 ℓ,α,τ for given 󰂃, chosen as described in [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF]. In comparison, a default implementation would evaluate the full cost function for each component of the gradient.

From decomposability to stochasticity

At each iteration k of the optimisation algorithm, the parameters are updated in an approximate fashion: θ (k) = θ (k-1) + ∆θ({∆r}, θ (k-1) ). The decomposition of f makes it possible minimize the cost function w.r.t. local parameters θ α,τ only at each iteration, and visit the parameters sequentially instead of jointly according to θ

(k) α,τ = θ (k-1) α,τ + ∆θ(∆r Bα,τ , θ (k-1)
Sα,τ ). Calculating the update ∆θ requires only subsets B α,τ and S α,τ of the full data and parameter sets, respectively.

If only gradient calculation and approximate line search are required, the smallest possible minibatch of data ∆r Bα,τ = {∆r α ′ ,τ ′ } (α ′ ,τ ′ )∈Bα,τ is such that: B α,τ = (R α ×{τ })∪({α}×T τ ), which is what we used in the main text. Here the notations {k} and {t} represent singletons and × the Cartesian (outer) product.

Similarly, the smallest possible subset of the parameters (subspace) that is accessed due to the regularizing priors is S α,τ = ∪ α ′ ,τ ′ ∈Bα,τ B α ′ ,τ ′ .

Altogether, we consider a minibatch sampling function:

B : k -→ (α, τ ) -→ {∆r α ′ ,τ ′ } (α ′ ,τ ′ )∈Bα,τ
This implies several possible optimization paths. If either the k -→ (α, τ ) or the (α, τ ) -→ {∆r α ′ ,τ ′ } (α ′ ,τ ′ )∈Bα,τ parts acts randomly, we thus define a stochastic approach.

Line search for descending the gradient in a subspace

As described in the Methods section, the update at iteration k takes the following form: ∆θ(∆r Bα,τ , θ

(k-1)
Sα,τ ) = s (k) p (k) , with descent direction given by p (k) = -H (k-1) ∇f α,τ (θ (k-1) ). The step size s (k) can be estimated in various ways. Here, in the context of minimizing f only w.r.t. the parameters θ α,τ (i.e. in the subspace D α,τ = Φ × {α} × {τ }, where Φ = {D, V }), s (k) is best estimated using an approximate line search. Indeed, for each cell (α, τ ), a line search operates in a subspace of small dimensionality, |D α,τ | (here two), and interestingly the same partial evaluation as for gradient calculation stands. This simplification comes from the fact that an approximate line search selects a step on basis of differences in f along the descent direction (Armijo's rule), i.e. with only the parameters in D that vary, and optionally of differences in ∇f (Wolfe's rule) in D as well.

Consequently, such an approximate line search is efficient and there is no need to determine s (k) using an adaptive rule, which is another departure of our algorithm from most of the other stochastic optimization procedures.

Inverse Hessian approximation

The inverse Hessian H was updated using the BFGS rule [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF]. Estimating the gradient and expressing the updates in small subspaces D . makes the updates to the inverse Hessian matrix similarly small (D . × D . ). This leads to the approximation for almost all parameter pairs that ∂f ∂θi∂θj (θ) ≈ 0. In our system this approximation is supported by the fact that separated regions in space and time have limited interactions regarding their parameters.

For example, the inference showcased in this article exhibits a block diagonal inverse Hessian matrix composed of 2×2 blocks in its stochastic version. We observed that the BFGS update rule for such a sparse inverse Hessian allowed faster convergence than a basic gradient descent approach, which corresponds to having the identity matrix replace H. This is illustrated Figure 4.

In complex cases, or when some domains contain limited numbers of points, H may not be definite positive. In this case we replace it by the identity matrix.

Correcting spatial oscillations in V

As Figure S4 in [13] shows, direct optimisation of the posterior distribution can lead to a potential energy landscape exhibiting spatial oscillations (of amplitude ≈ 0.15 k B T for a diffusivity of D ≈ 0.1 µm 2 s -1 ) even in the absence of forces.

The posterior probability [Eq. (9)] only depends on the potential energy V through its spatial gradient ∇V . A gradient is insensitive to spatial oscillations on a regular mesh at the Nyquist frequency -∇V is always evaluated on top of a hill or at the bottom of a trough and the variation in-between is just not sampled. This undesirable phenomenon, limited in full optimization, was amplified when locally and sequentially updating the parameters since the approximate line search generates noisy updates that are even noisier when the step sizes are driven by fewer parameters.

Regularising priors lessen the problem as they add a damping factor. However, when attempting to detect subtle effects and features in the map, it is not desirable to increase regularisation beyond what is needed to attenuate statistical errors. We instead introduced two modifications of the approximate line search, which also make the algorithm more robust in general:

• setting an upper bound on individual updates,

• adapting Wolfe's rule in the approximate line search to prevent gradient reversal.

An upper bound on the infinite-norm of the update allows to scale the initial update considered in the line search (for our stochastic algorithm, the infinite norm is equal to the maximum of the absolute values of D α,τ and V α,τ ). This mostly makes the search faster and better uses a limited number of iterations (default is 10), but can also prevent large updates. We set this upper bound equal to 2k B T for V and 2µm 2 s -1 for D based on physical considerations, though the precise value of the bound is not important as its main purpose is to avoid diverging updates.

The standard strong Wolfe's rule [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] aims at making the norm of the projected gradient decrease (with iteration index k shown as subscript for convenience):

cp T k ∇f (θ k-1 ) ≤ -p T k ∇f (θ k-1 +s k p k ) ≤ -cp T k ∇f (θ k-1 ) , (17)
with c = 0.9. We propose instead:

cp T k ∇f (θ k-1 ) ≤ -p T k ∇f (θ k-1 + s k p k ) , (18)
with c = 0.5 (or lower).

Our first motivation here is to bound more tightly the projected gradient p T k ∇f (θ k-1 + s k p k ) when it points opposite to p T k ∇f (θ k-1 ). Indeed, a gradient reversal may be the signature of a large step and cause oscillations around the optimal parameter value, especially with local updates. With c = 0.5, a reversed gradient is admitted only if it is half the size of the initial gradient, along the descent direction.

A second motivation in dropping the weak part of Wolfe's rule deals with non-monotonic properties in the descent subspaces. Empirically, situations revealed cases where the projected gradient was only increasing along the descent direction, while the candidate updates complied with the Armijo rule [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] (decrease in f ). As long as the Armijo rule is ensured, if the gradient keeps on pointing in the same direction, we only need to carry out more steps. Especially, if the projected gradient increases in absolute value, then a small step is preferable in order to stay in smooth conditions, so that the inverse Hessian matrix can be properly updated.

Parallelisation of the optimisation

The main constraint in parallelising the computation performed on minibatches is to ensure non-conflicting access to resources during execution of the various tasks. We implemented the parallel procedure the following way.

A scheduler process operates the k -→ (α, τ ) sampling function and dispatches single iterations (tasks) to a pool of workers. Each worker will process an iteration and consequently update the parameters of a given spatio-temporal domain, while another worker will simultaneously process another iteration. Parameters involved in the computations run by a worker are not modified during computation by any other workers.

The parameters involved at iteration k are identified by the following function: S : k -→ (α, τ ) -→ θ Sα,τ , where θ Sα,τ = {θ α ′ ,τ ′ } (α ′ ,τ ′ )∈Sα,τ . Each time the scheduler assigns an iteration to a domain (α, τ ), it can lock the corresponding S α,τ parameters. The later iterations that require access to locked parameters will be postponed until the locked parameters are released. Finally, note that to allow for more parallelisation on small maps, we locked B α,τ instead of S α,τ .

More precisely, each worker maintains its own copy of the parameter vector θ. The workers synchronise together, exchanging local parameter updates θ

(k) = f 󰀃 θ (k-1) 󰀄 -f 󰀃 θ (k) 󰀄 .
Convergence is stated once ∆f drops below some small tolerance for 90% of the iterations in an epoch of ΩT iterations, where ΩT is the number of domains in the map.

To ensure that all domains are updated, they are each sampled once, in random order, during each epoch. To avoid potential latencies between epochs, the next epoch begins as soon as all domains have been sampled; the scheduler does not wait for the workers to complete the corresponding updates.

TIME REGULARISATION AND INITIAL

VALUES FOR V

In the context of inference of a potential energy field, time regularization is especially useful in making potential energies comparable between different time segments. Nothing in the underlying physics fixes the average value of the effective potential, creating complexities when evaluating the effective potential at different time points.

In InferenceMAP [24], the potential fields were initialized to match the corresponding equilibrium distribution to the recorded point density. When temporal regularisation is imposed, direct optimisation does not necessarily constrain the effective potential (in non interacting space domains) for small values of λ t . To make temporal smoothing less sensitive to the value of λ t and to ensure similar levels between time segments, we initialised V to 0 instead, for inferences with temporal regularisation.

FIG. 2 .

 2 FIG.2. Inferred potential energy and diffusivity landscapes in the static, low-dimensional example (558 parameters in total).

FIG. 6 .

 6 FIG. 6. Effects of temporal regularisation on inferred maps for an example consisting in two static and one temporally varying potential wells. The centers of the wells are all placed on the central horizontal line of the maps, with the temporally varying well in the center and the static wells to each side (see top rows). (Top rows) Two examples of inferred spatio-temporal maps for different strengths of the temporal regularisation (log 10 λt = -2 and log 10 λt = +1, respectively). (Bottom row) Temporal (left) and spatial (right) profiles of the inferred potential V for a selection of different values of the hyperparameter λt controlling temporal regularisation (logarithmic scale). (Bottom, left) The temporal profile shows the value in the central space domain (see top rows) as function of time. (Bottom, right) The spatial profile shows the values in the spatial domains lying on the horizontal line going through the central domain at the time when the depth of the potential well is maximal.

FIG. 7 .

 7 FIG. 7. Potential energy maps of Gag proteins in part of a T cell (127 578 parameters in total and 126 022 recorded localisations). The second row shows the full inferred map. The rows above and below show magnifications of spatial regions of interest. White domains were not visited by any molecule and not represented (no associated parameter).

  of the inverse Hessian matrix. All workers send status information to the scheduler, this merely includes the decrease in f : ∆f

SUPPLEMENTARY INFORMATION SAMPLING THE POSTERIOR DISTRIBUTION

Here, we relied on stochastic optimisation to access the maximum a posteriori (MAP) values of the map parameters. Sampling of the posterior distribution around the MAP can be performed using a similar procedure, i.e. stochastic gradient Langevin dynamics [START_REF] Welling | Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML'11[END_REF]. It uses the following update rule:


Lorenzo Guglielmetti

Nicolas Veziris

Alexandra Aubry

Florence Brossier

Christine Bernard

Wladimir Sougakoff

Vincent Jarlier

MD Jérôme Robert
email: jerome.robert@aphp.fr

Risk-factors for extensively drug resistance in multidrug-resistant tuberculosis cases: a case-case study AUTHORS

Settings. Extensively drug-resistant (XDR) tuberculosis (TB) identification may be delayed because of the lack of availability of molecular testing to second-line drugs. Hence, early suspicion of XDR-TB is necessary to avoid developing further drug resistance.

Objective. We sought to identify the characteristics associated with XDR-TB among multidrug-resistant (MDR) TB cases, before the availability of second-line drug susceptibility testing (DST) results.

Methods.

All MDR-TB cases with available second-line DST results recorded in France from 1998 to 2013 were classified as simple MDR (no resistance to fluoroquinolones and secondline injectables), pre-XDR (resistance to fluoroquinolones or second-line injectables), and XDR cases (resistance to both).

Results.

A total of 833 MDR-TB cases were analysed, including 168 (20%) pre-XDR and 62 (7%) XDR-TB cases. A prior history of treatment was acknowledged for 41% of the cases, and 12% were HIV-positive. Characteristics independently associated with XDR-TB were foreign birth (odds ratio: 9.5), prior TB treatment (OR: 2.6), smear positivity (OR: 4.5), and ethambutol resistance (OR: 9.1). Characteristics independently associated with pre-XDR-TB when compared to simple MDR-TB cases were male gender (OR: 1.6), birth in Europe (OR: 2.6), and ethambutol resistance (OR: 1.9).

Conclusion

. The presence of clinical or bacteriological characteristics associated with XDR-TB should lead to rapid molecular testing for resistance to second-line drugs before starting a tailored treatment. (OR : 1,9). Conclusion : Chez un malade porteur de TB MDR, la présence de caractéristiques cliniques ou bactériologiques associées à la TB XDR doit faire réaliser des S2L sans délai avant de débuter un traitement adapté individuellement.

Résumé

Cadre : L'identification des cas de tuberculose (TB) à bacilles ultra-résistants (XDR) aux antibiotiques peut être retardée par manque d'accès aux tests de sensibilité aux antituberculeux de seconde ligne (S2L). Afin de ne pas risquer d'aggraver l'étendue de la résistance, il est important d'identifier ces cas XDR sur d'autres critères.

Objectif : Identifier les caractéristiques associées à la TB XDR parmi les cas identifiés de TB à bacilles multirésistants (MDR), avant disponibilité des résultats des S2L.

Méthodes : Tous les cas de TB MDR avec des résultats de S2L diagnostiqués en France de 1998 à 2013 ont été classés en cas MDR simple (pas de résistance aux fluoroquinolones ni aux antibiotiques injectables de seconde ligne), pré-XDR (résistance à un seul de ces deux groupe d'antibiotiques), et XDR (résistance aux deux groupes d'antibiotiques).

Résultats : AU total, 833 cas de TB MDR ont été inclus dans l'analyse, dont 168 (20%) cas pré-XDR et 62 (7%) cas XDR. Des antécédents de traitement antituberculeux ont été identifiés chez 41% du total des cas et 12% étaient séropositifs pour le VIH. Les caractéristiques significativement associées à la TB XDR étaient la naissance à l'étranger (odds ratio : 9,5), les antécédents de traitements (OR : 2,6), un examen microscopique positif (OR : 4,5) et la résistance à l'éthambutol (OR : 9,1). Par comparaison au cas MDR simples, les caractéristiques significativement associées à la TB pré-XDR étaient le sexe masculin (OR : 1,6), la naissance en Europe hors de France (OR : 2,6) et la résistance à l'éthambutol

INTRODUCTION

France is classified as a low incidence country for tuberculosis (TB) [START_REF] Who | Global tuberculosis report[END_REF]. Each year approximately 6000 cases of tuberculosis are diagnosed, representing a prevalence rate of 8-9 per 100 000, and between 50 to 100 cases or 1% to 2% of the total are multidrug-resistant (MDR) (i.e., resistant to isoniazid and rifampicin). A striking increase in the number of MDR-TB cases due to migration from Eastern Europe and the Caucasus has been recently reported in France (2,3), while the MDR-TB prevalence had been stable during the previous years (4).

This changing epidemiology has been observed in other low-incidence countries in Western

Europe, where the phenomenon poses a substantial threat to the World Health Organization (WHO) goal of TB elimination [START_REF] Falzon | Drug resistance and tuberculosis elimination in low-incidence countries[END_REF]. It has been estimated that MDR-TB/rifampicin-resistant cases occurring in the WHO European region, mostly in its Eastern part, account for more than one fifth of the total MDR-TB/rifampicin-resistant cases estimated in the world in 2015 [START_REF] Who | Global tuberculosis report[END_REF].

Among MDR-TB cases recently reported in western European countries, increasing numbers fulfil the definition of extensively drug-resistant (XDR) TB, i.e. MDR isolates with additional resistance to any fluoroquinolone and any second-line injectable drug (SLI). Despite prolonged treatment regimens, which according to WHO recommendations can exceed 20 months (6), treatment success was achieved in only 64% of MDR-TB cases in a large metaanalysis [START_REF] Falzon | Resistance to fluoroquinolones and second-line injectable drugs: impact on multidrug-resistant TB outcomes[END_REF]. The proportion of successful outcomes was lower for pre-XDR TB, defined as MDR-TB with additional resistance to any fluoroquinolone or any SLI, and further decreased to 40% for XDR-TB cases [START_REF] Falzon | Resistance to fluoroquinolones and second-line injectable drugs: impact on multidrug-resistant TB outcomes[END_REF]. Moreover, the prognosis of XDR-TB cases has been shown to be worse for patients harbouring isolates with additional resistance to other first-and secondline drugs [START_REF] Migliori | Drug resistance beyond extensively drug-resistant tuberculosis: individual patient data meta-analysis[END_REF]. The mortality associated with XDR-TB varies according to the degree of HIV prevalence of the population, while in all settings a remarkable number of these subjects are discharged in the community, which contributes to the spread of the disease [START_REF] Pietersen | Long-term outcomes of patients with extensively drug-resistant tuberculosis in South Africa: a cohort study[END_REF][START_REF] Keshavjee | Treatment of extensively drug-resistant tuberculosis in Tomsk, Russia: a retrospective cohort study[END_REF].

The increasing burden of XDR-TB, together with its difficult management and unsatisfactory outcome, emphasizes the importance of timely diagnosis and treatment initiation for XDR-TB patients. However, drug susceptibility tests (DST) are not available in all settings worldwide [START_REF] Who | Global tuberculosis report[END_REF]. The development of automated molecular assays for DST to first line drugs may facilitate the diagnosis of MDR tuberculosis but is currently not implemented for all TB cases. Rapid molecular testing for second-line drugs is still restricted to selected laboratories, even in highincome countries. It is therefore of interest to identify the characteristics associated with additional resistance in MDR-TB patients for early recognition and better early management in order to avoid increasing the initial burden of resistance. The major objective of our study is to identify epidemiological and bacteriological factors associated with XDR-TB when compared to non-XDR MDR-TB patients diagnosed in France.

METHODS

MDR-TB surveillance in France

Since 1992, the national Reference Centre for Mycobacteria (NRC) pilots a network of 250 laboratories, which includes all centres performing mycobacterial culture in France. Among them, approximately a third are performing drug susceptibility tests (DST) for first line drugs by phenotypic methods and less than 15 (6%) are using molecular biology to test for resistance to second-line drugs. Each year, the laboratories report the number of MDR-TB cases to the NRC, which confirm the MDR status and retrospectively collects additional information on MDR-TB patients by sending a standardized data collection form. Initial clinical and demographical data are collected for each MDR-TB patient, including country of birth, HIV status, previous TB treatment history, and microbiology. Finally, MDR Mycobacterium tuberculosis complex isolates are sent to the NRC laboratory in order to perform comprehensive DST for first and second-line anti-TB drugs. From 1992 to 2004, the second-line DST were performed at either Pitié-Salpêtrière Hospital or the Pasteur Institute, in Paris. Since 2005, this task was entirely taken over by the Pitié-Salpêtrière NRC.

Phenotypic DST for first-and second-line drugs were performed by using the proportion method on Löwenstein-Jensen medium according to WHO recommendations [START_REF] Who | Guidelines for Surveillance of Drug Resistance in tuberculosis[END_REF].

Case-case studies

Since DST for fluoroquinolones and SLI started to be routinely performed in 1998, only MDR-TB cases reported from 1998 to 2013 were included in the analysis. Patients reported during multiple years, i.e. culture-positive chronic cases, were considered only at the first occurrence. Duplicates were eliminated before data anonymization. Patients with missing DST result for fluoroquinolones and SLI were excluded from the analysis.

In order to compare characteristics of MDR-TB patients harbouring different resistance patterns to fluoroquinolones and SLI, we performed different nested case-case studies within the cohort. In the main analysis, patients harbouring XDR isolates were compared to all the other cases harbouring non-XDR MDR isolates. Two other case-case analyses were conducted, comparing i) cases with pre-XDR isolates to those harbouring MDR isolates with no additional resistance to any fluoroquinolone and any SLI (defined as simple-MDR), and ii) XDR to pre-XDR patients.

Data analysis

Medians with interquartile ranges (IQR) were calculated for continuous variables, and frequency distributions were tabulated for categorical variables. Groups were compared using the two-sample Wilcoxon-Mann-Whitney test for continuous data and Fisher's exact test for categorical data. In order to assess the independent association of predictive variables, logistic regression modelling was performed. For each case-case analysis, variables associated to a group of patients under a cut-off of p-value <0.2 in univariate analysis were added to a multivariable model. The variables of resistance to ethambutol and streptomycin were considered for multivariable analysis, as they could be of interest to predict additional resistance in MDR-TB cases. Pyrazinamide DST status was not included because of the difficulty to interpret its result and consequently a high number of missing data. Patients with unknown HIV status were considered as HIV-negative and patients with unknown place of birth were classified as foreign-born. With regard to the other variables, patients with missing data were excluded from the analysis. Backward selection procedure was performed, with age and sex included by default in all models. The foreign-born status or WHO regions of birth were tested alternatively in the models. Birth in France was considered as the reference in all models. The best-fitted model was chosen on the basis of epidemiological and statistical criteria, such as the Hosmer-Lemeshow goodness-of-fit test and the Akaike Information Criterion. Interaction terms were included in the models and dropped if non-significant. Only the model with the least number of variables is reported for each case-case analysis. Odds ratios (OR) were calculated with their 95% confidence intervals (CI). P-values < 0.05 were considered as statistically significant. STATA, version 11.1 (StataCorp, College Station, Texas, USA) was used for statistical analysis.

Ethics

Data were routinely collected during standard care and abstracted from medical files. No personal identifiers were recorded. Therefore, according to French regulations, the study did not require institutional review board or patients' approvals.

RESULTS

Overall, 983 MDR-TB cases have been recorded in France during the 16-year study period, ranging from 31 to 91 per year. Among them, 93 (9.5%) cases were discarded because of duplicate declarations and 57 (5.8%) because DST results were not available (Figure 1). The baseline characteristics of the 833 MDR-TB patients finally included in the study are shown in Table 1. Patients were mostly male (64%), born in Sub-Saharan Africa (30%) or in European countries excluding France (28%). Among the latter, a large majority of patients (91%) were from the 18 high-priority countries (HPC) identified by WHO Regional Office for Europe's Stop TB Strategy.

French-born patients accounted for only 15% of the cases, and the country of birth was unknown for 23 (2.8%) patients. A majority of patients had pulmonary TB (90%), and 56% had positive sputum-smear examination. A total of 41% had a prior history for TB treatment.

A total of 12% were HIV-positive, 76% HIV-negative, and 12% had unknown HIV status.

The median age at diagnosis was 32 years (IQR: 25-43).

A total of 603 (72%) patients had isolates susceptible to both fluoroquinolones and SLI, 168 (20%) pre-XDR isolates, and 62 (7%) XDR isolates. Additional resistance was frequent for streptomycin (69%) and ethambutol (47%).

Comparison of XDR-TB cases to all other MDR-TB cases

When comparing the 62 XDR-TB to the 771 non-XDR MDR-TB cases, XDR-TB was associated with an age at diagnosis of 35 to 44 years (Table 1). XDR-TB patients were significantly more likely to be born outside France and particularly in Europe (OR: 35.6; 95% CI 13.9-91.6), to have a prior history of TB treatment (OR: 3.7; 2.0-6.9), and to be sputum smear-positive (OR: 3.5; 1.8-6.9). No significant association was observed with sex, disease localisation and HIV status. XDR isolates were significantly more likely to harbour resistance to streptomycin and ethambutol (Table 1).

After multivariable logistic regression analysis, foreign birth (OR 9.5; CI: 1.2-75.9), sputum smear positivity (OR 4.5; CI: 2.0-10.1), prior TB treatment (OR 2.6; CI: 1.4-5.0), and resistance to ethambutol (OR 9.1; CI: 4.0-20.7) remained independently associated with the risk of XDR-TB (Tables 2, model 1). When replacing foreign birth by regions of birth in a second model (Table 2, model 2), birth in the WHO European region (OR 21.0; CI: 7.0-63.1) was significantly associated with XDR-TB, while no other region of birth remained associated with XDR-TB status.

Comparison of pre-XDR TB cases to simple-MDR-TB cases

In univariate analysis, pre-XDR TB cases were significantly more likely than simple-MDR TB cases to be male, born in the European region, born in Africa excluding the Sub-Saharan region, to have prior TB treatment and to harbour isolates with additional resistance to ethambutol and streptomycin (data not shown). In multivariable analysis, characteristics significantly associated with pre-XDR TB as compared to simple-MDR TB were male gender (OR 1.6; CI: 1.0-2.6), birth in Europe (OR 2.6; CI: 1.3-5.1), and ethambutol resistance (OR 1.9; CI: 1.4-2.7) (Table 3, model 1).

Comparison of pre-XDR TB cases to XDR-TB cases

In univariate analysis, XDR-TB cases were significantly more likely than pre-XDR TB cases to be born in the European region, or in Africa excluding the Sub-Saharan region, to have a positive sputum smear examination, to have prior TB treatment, and to harbour isolates with additional resistance to ethambutol (data not shown). In multivariable analysis, characteristics remaining significantly associated with XDR-TB as compared to pre-XDR TB were birth in Europe (OR 8.1; CI: 3.1-21.3), sputum smear positivity (OR 5.2; CI: 2.9-9.3), and ethambutol resistance (OR 3.9; CI: 1.7-8.8) (Table 3, model 2).

DISCUSSION

We compared the characteristics of patients diagnosed with MDR-TB or XDR-TB over 16 years in France. The results of our study delineate a few characteristics associated with resistance beyond MDR.

First, foreign birth was associated with XDR-TB, and with higher resistance in both intermediate steps between MDR-TB and XDR-TB. In low TB incidence countries, the predominance of foreign-born patients among MDR-TB cases is a common finding [START_REF] Eker | Multidrug-and Extensively Drug-Resistant Tuberculosis, Germany[END_REF][START_REF] Schecter | Linezolid in the Treatment of Multidrug-Resistant Tuberculosis[END_REF][START_REF] Anderson | Treatment outcome of multidrug resistant tuberculosis in the United Kingdom: retrospective-prospective cohort study from 2004 to 2007[END_REF][START_REF] Van Altena | Highly successful treatment outcome of multidrug-resistant tuberculosis in the Netherlands, 2000-2009[END_REF].

In our study, the highest risk group was those of patients born in countries with high incidence of MDR/XDR-TB, such as European HPC. This phenomenon has been previously described in a country-wide cohort of TB patients in Canada [START_REF] Minion | Multidrug and Extensively Drug-resistant Tuberculosis in Canada 1997-2008: Demographic and Disease Characteristics[END_REF]. Second, our results have shown that previous TB treatment was significantly associated with XDR-TB. This is a wellknown risk factor, which has already been described in different settings in association with both MDR-TB [START_REF] Andrews | Predictors of Multidrugand Extensively Drug-Resistant Tuberculosis in a High HIV Prevalence Community[END_REF][START_REF] Vashakidze | Prevalence and risk factors for drug resistance among hospitalized tuberculosis patients in Georgia[END_REF][START_REF] Skrahina | Alarming levels of drug-resistant tuberculosis in Belarus: results of a survey in Minsk[END_REF][START_REF] Liang | Factors contributing to the high prevalence of multidrug-resistant tuberculosis: a study from China[END_REF] and XDR-TB [START_REF] Shin | Development of Extensively Drug-resistant Tuberculosis during Multidrug-resistant Tuberculosis Treatment[END_REF][START_REF] Kliiman | Predictors of Extensively Drug-Resistant Pulmonary Tuberculosis[END_REF][START_REF] Jeon | Extensively Drug-Resistant Tuberculosis in South Korea: Risk Factors and Treatment Outcomes among Patients at a Tertiary Referral Hospital[END_REF][START_REF] Dalton | Prevalence of and risk factors for resistance to second-line drugs in people with multidrug-resistant tuberculosis in eight countries: a prospective cohort study[END_REF]. Third, smear positivity was also independently associated with XDR-TB in our cohort, in the comparison with both MDR-TB and pre-XDR TB strains. This could be linked to the fact that most XDR-TB patients in our cohort had a long previous treatment histories and advanced pulmonary lesions. This finding, however, is not confirmed by other similar studies [START_REF] Dalton | Prevalence of and risk factors for resistance to second-line drugs in people with multidrug-resistant tuberculosis in eight countries: a prospective cohort study[END_REF][START_REF] Balabanova | Risk factors for drug-resistant tuberculosis patients in Lithuania, 2002-2008[END_REF]. This may also reflect high bacterial load, which is directly linked to a higher likelihood of resistant mutants in the lesions. Fourth, XDR-TB cases were associated with significantly higher rates of ethambutol resistance, and the same finding was reported for more resistant strains in both intermediate analyses. Higher rates of ethambutol resistance have been consistently reported in XDR-TB patients [START_REF] Keshavjee | Treatment of extensively drug-resistant tuberculosis in Tomsk, Russia: a retrospective cohort study[END_REF][START_REF] Mitnick | Comprehensive treatment of extensively drug-resistant tuberculosis[END_REF][START_REF] Günther | Multidrug-Resistant Tuberculosis in Europe, 2010-2011[END_REF], and may have been a facilitating factor for the progressive selection of resistance to second-line drugs. Hence, rapid availability of ethambutol DST could help in the management of MDR-TB patients. However, the performance of the first version of the most commonly used line probe assay for the detection of resistance to second-line drugs was suboptimal for ethambutol, probably due to the lack of testing in regions outside embB codon 306 [START_REF] Brossier | Molecular Analysis of the embCAB Locus and embR Gene Involved in Ethambutol Resistance in Clinical Isolates of Mycobacterium tuberculosis in France[END_REF]. Unfortunately, the second version of this test does not include ethambutol testing at all.

In addition, male sex was linked to pre-XDR TB when compared to MDR-TB, as previously reported in a study in the Russian Federation [START_REF] Shin | Development of Extensively Drug-resistant Tuberculosis during Multidrug-resistant Tuberculosis Treatment[END_REF], but in contrast to a multicentre study where female sex was linked to XDR-TB [START_REF] Dalton | Prevalence of and risk factors for resistance to second-line drugs in people with multidrug-resistant tuberculosis in eight countries: a prospective cohort study[END_REF]. Finally, we did not find a significant association between XDR-TB and HIV infection, in concordance with a large study including countries with variable HIV burden [START_REF] Dalton | Prevalence of and risk factors for resistance to second-line drugs in people with multidrug-resistant tuberculosis in eight countries: a prospective cohort study[END_REF]. In contrast, a link between HIV and XDR-TB has been shown in a high-incidence country like South Africa [START_REF] Andrews | Predictors of Multidrugand Extensively Drug-Resistant Tuberculosis in a High HIV Prevalence Community[END_REF].

Our study has some limitations. Risk factors that have been consistently described in literature, such as alcohol and drug abuse, history of imprisonment, and previous treatment with second-line drugs, were not captured in the routine data collection form and therefore not analysed in our study. In addition, the epidemiological changes that occurred in the MDR-TB population in France over time [START_REF] Bernard | A surge of MDR and XDR tuberculosis in France among patients born in the Former Soviet Union[END_REF] could have influenced our findings, and only a small proportion (7%) of patients included in the study were XDR-TB. However, both these secular trends and the low rates of XDR-TB cases are common to most low-incidence countries [START_REF]Tuberculosis surveillance and monitoring in Europe 2016[END_REF], and therefore do not affect the generalizability of our results to these settings. Finally, a limited number of patients had missing data for the key variables (HIV status and country of birth), which were managed by making assumptions.

In conclusion, we identified factors associated with XDR-TB among patients diagnosed with MDR-TB in a low-incidence setting. These characteristics should increase the level of suspicion for XDR-TB among MDR-TB cases and warrant testing with rapid molecular assays for second-line drugs. Indeed, a new version of a line probe assay has been shown to be improved [START_REF] Brossier | Performance of the New Version (v2.0) of the GenoType MTBDR sl Test for Detection of Resistance to Second-Line Drugs in Multidrug-Resistant Mycobacterium tuberculosis Complex Strains. Land GA[END_REF] and has recently been endorsed by the WHO [START_REF] Who | (SLIDs): policy guidance[END_REF]. The combination of implementing rapid diagnostic tools and identifying high-risk groups may lead to promptly detect XDR-TB cases and to design effective empirical regimens while waiting for comprehensive DST results. Finally, new shorter MDR-TB treatment regimens should be avoided in patients harbouring such characteristics.

Figure 1 .

 1 Figure 1. Flow-chart of the study 304

Table 1 .

 1 Characteristics of 833 multidrug-resistant tuberculosis (MDR-TB) cases according to extensively drug-resistant (XDR-TB) status.

		All cases	XDR-TB isolates	OR	p
	Categorical variables	(n=833)	Yes	No	(CI 95%)
			(n=62)	(n=771)
		n (%)	n (%)	n (%)
	Sex Male	533 (64)	45 (74)	488 (64)	1.6 (0.8-2.9)	0.16
	Age*:
	≤14	26 (3)	1 (2)	25 (3)	1
	15-24	167 (20)	9 (14)	158 (20)	1.4 (0.5-4.0)	0.51
	25-34	298 (36)	19 (31)	279 (36)	1.7 (0.7-4.2)	0.25
	35-44	165 (20)	21 (34)	144 (19)	3.6 (1.5-8.7) <0.01
	45-64	133 (16)	10 (16)	123 (16)	2.0 (0.6-6.4)	0.23
	≥65	38 (5)	2 (3)	36 (5)	1.4 (0.3-5.4)	0.64
	Foreign-birth or unknown country*	710 (85)	61 (98)	649 (84) 11.5 (4.2-31.6) 0.001
	Region of birth*:
	France	123 (15)	1 (2)	122 (16)	1	-
	Europe (France excluded)	230 (28)	52 (84)	178 (24) 35.6 (13.9-91.6) <0.001
	North Africa	97 (12)	4 (6)	93 (12)	5.2 (0.9-30.4) 0.06
	Sub-Saharan Africa	252 (30)	3 (5)	249 (33)	1.5 (0.3-6.3)	0.61
	Others	108 (13)	2 (3)	106 (14)	2.3 (0.6-9.2)	0.24
	Pulmonary tuberculosis	747 (90)	60 (97)	687 (90)	3.3 (1.1-10.1)	0.11
	Smear positive for pulmonary cases	464 (56)	51 (82)	413 (57)	3.5 (1.8-6.9) <0.001
	HIV-positive*	101 (12)	4 (6)	97 (13)	0.5 (0.1-2.0)	0.22
	Prior tuberculosis treatment	343 (41)	43 (73)	300 (42)	3.7 (2.0-6.9) <0.001
	Resistance to:
	Ethambutol	389 (47)	55 (89)	334 (44) 10.2 (3.6-28.9) <0.001
	Streptomycin	577 (69)	53 (85)	524 (69)	2.6 (1.3-5.3)	0.006
	OR, odds ratio; CI, confidence interval
	*6 patients with missing age,23 with unknown region of birth, and 97 with unknown HIV status

Table 2 .

 2 Multivariable logistic regression analysis assessing characteristics associated with extensively drug-resistant tuberculosis in a cohort of multidrug-resistant tuberculosis cases. The country of birth of patients was entered into the model after categorization according to foreign-born versus French-born (Model 1) or according to the birth in the European region versus birth outside the European region and French-born (Model 2). The European region was defined according to the World Health Organisation; France was excluded from this group.

		Model 1		Model 2
	Categorical variables	OR	p	OR	p
		(CI 95%)		(CI 95%)
	Foreign birth	9.5 (1.2 -75.9)	0.034	-
	Born in the European region (France excluded)	-		21.0 (7.0 -63.1) <0.001
	Sputum smear-positivity	4.5 (2.0 -10.1)	<0.001	3.8 (2.1 -6.8)	<0.001
	Prior tuberculosis treatment	2.6 (1.4 -5.0)	0.004	2.2 (1.2 -3.9)	0.018
	Additional resistance to ethambutol	9.1 (4.0 -20.7)	<0.001	6.1 (2.3 -16.4)	<0.001
	OR, odds ratio; CI, confidence interval.

Table 3 .

 3 Multivariable logistic regression analysis assessing characteristics associated with tuberculosis (TB) due to multidrug-resistant (MDR) isolates with additional resistance to any fluoroquinolone or any second-line injectable drug (pre-XDR) as compared to MDR-TB with no additional resistance to any fluoroquinolone or second-line injectable drug (Model 1) and characteristics associated with extensively drug-resistant TB as compared to Pre-XDR TB (Model 2).The European region was defined according to the World Health Organisation; France was excluded from this group and included as itself in the models.

		Model 1		Model 2
		MDR-TB vs. pre-XDR TB	pre-XDR TB vs. XDR-TB
	Categorical variables	OR	p	OR	p
		(IC 95%)		(IC 95%)
	Sex, male	1.6 (1.0 -2.6)	0.033	-
	Born in the European region (France excluded)	2.6 (1.3 -5.1)	0.006	8.1 (3.1 -21.3)	<0.001
	Sputum smear-positive pulmonary tuberculosis	-		5.2 (2.9 -9.3)	<0.001
	Strain with additional resistance to Ethambutol	1.9 (1.4 -2.7)	<0.001	3.9 (1.7 -8.8)	<0.001
	CI, confidence interval; OR, odds ratio.

ACKNOWLEDGEMENTS

The authors thank all microbiologists that participated in the NRC Laboratory Network for MDR-TB surveillance and Cristina Gutierrez for data checking during the 1998-2004 period.

We are also thankful to the clinicians who shared information.

JR conceived and designed the study

LG & JR performed the statistical analysis, analysed the results and drafted the manuscript

LG, NV, AA, FB, CB, WS, VJ & JR collected the data, participated in the interpretation of the results, and critically revised the manuscript.

NV, AA, FB, CB & WS were in charge of drug susceptibility tests.

All authors gave final approval of the current version to be published BIBLIOGRAPHY