Michal Valko
email: michal@cs.pitt.edu

Milos Hauskrecht

Distance Metric Learning for Conditional Anomaly Detection

Anomaly detection methods can be very useful in identifying unusual or interesting patterns in data. A recently proposed conditional anomaly detection framework extends anomaly detection to the problem of identifying anomalous patterns on a subset of attributes in the data. The anomaly always depends (is conditioned) on the value of remaining attributes. The work presented in this paper focuses on instance-based methods for detecting conditional anomalies. The methods depend heavily on the distance metric that lets us identify examples in the dataset that are most critical for detecting the anomaly. To optimize the performance of the anomaly detection methods we explore and study metric learning methods. We evaluate the quality of our methods on the Pneumonia PORT dataset by detecting unusual admission decisions for patients with the community-acquired pneumonia. The results of our metric learning methods show an improved detection performance over standard distance metrics, which is very promising for building automated anomaly detection systems for variety of intelligent monitoring applications.

Introduction

Anomaly detection methods can be very useful in identifying interesting or concerning events. Typical anomaly detection attempts to identify unusual data instances that deviate from the majority of examples in the dataset. Such instances indicate anomalous (out of ordinary) circumstances, for example, a network attack [START_REF] Eskin | Anomaly detection over noisy data using learned probability distributions[END_REF] or a disease outbreak [START_REF] Wong | Bayesian network anomaly pattern detection for disease outbreaks[END_REF]. In this work, we study conditional anomaly detection framework that extends standard anomaly detection by identifying partial patterns in data instances that are anomalous with respect to the remaining data features. Such a framework has been successfully applied to identify unusual patient-management decisions made for patients suffering from different conditions [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF].

Data attributes (features) in the conditional anomaly detection are divided into two disjoint groups: context (or condition) attributes C and target attributes A. The conditional Copyright c 2011, Association for the Advancement of Artificial Intelligence (www.aaai.org). All rights reserved.

anomaly methods then attempt to identify anomalies in target attributes A with respect to context C. The conditional aspect allows us to identify patterns that are typical in one context but anomalous in the other. To illustrate the potential of the method assume two patients with different conditions are given the same drug. In one of these conditions the administration of the drug is normal, but for the other one it is unusual. The conditional anomaly detection methods with target variables corresponding to the treatment should be able to identify the anomaly in the treatment.

The conditional anomaly detection method evaluates and identifies anomalies one data example at the time. To make an anomaly call for a data instance x, [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF] proposed a probabilistic predictive model M that aims to capture stochastic dependencies among the target and context attributes. The predictive model defines a conditional probability distribution p(A|C) of target attributes given the values of context variables. Given the predictive model, the anomaly call for a data instance x is made if the probability of the target attributes observed in x is small.

A predictive probabilistic model used for detection purposes can be built in different ways. In this paper, we focus on instance-based approaches. The instance-based methods do not try to learn a universal predictive model for all possible instances at the same time, instead the model is optimized for every data instance x individually. The instancespecific model M x may provide a better option if the predictive model is less complex and the dataset is small.

Instance-specific models often rely on a distance metric that aims to pick examples most relevant for the prediction of x. However, the question of what is the best distance metric to reflect the relevancy of the example to the prediction is the most challenging part of the task. Standard metrics such as Euclidean or Mahalanobis metrics are not the best for the anomaly detection task since they may be biased by feature duplicates or features that are irrelevant for predicting target attributes. Thus, instead of choosing one of the standard distance metrics we investigate and test metric-learning methods that let us adapt predictive models to specifics of the currently evaluated example x.

We explore two metric-learning methods that were originally used for building non-parametric classification models. The first method is NCA [START_REF] Goldberger | Neighbourhood components analysis[END_REF]). The method adjusts the parameters of the generalized distance metric so that the accuracy of the associated nearest neighbor classifier is optimized. The second method, RCA [START_REF] Bar-Hillel | Learning a mahalanobis metric from equivalence constraints[END_REF], optimizes mutual information between the distribution in the original and the transformed space with restriction that distances between same class cases do not exceed a fixed threshold.

To evaluate the quality of metric learning methods in anomaly detection we apply them to the problem of identification of unusual patient-management decisions, more specifically, to the problem of detection of unusual hospitalization patterns for patients with the community acquired pneumonia. We show that on this problem metric learning approaches outperform standard distance metrics.

Methodology Conditional anomaly detection

In anomaly detection, we are interested in detecting an unusual data pattern the occurrence of which deviates from patterns seen for other examples. In the conditional anomaly, a partial data pattern is evaluated in context of other data variables and their values. Briefly, the data attributes (features) are divided into two disjoint groups: context (or condition) attributes C and the target attributes A. The objective of conditional anomaly detection methods is to identify anomalies in target attributes A with respect to context attributes C.

Let E = {x 1 , x 2 , . . . , x n } be a set of examples in the dataset and let x be an example we want to analyze and determine if it is conditionally anomalous with respect to examples in E. The context of the example x is defined by the projection of x to context attributes C, which we denote by C(x). Similarly, A(x) denotes the projection of x to target attributes.

Our goal is to identify the anomaly in x with respect to examples in the dataset E. The examples and their relation to x can be captured indirectly by an auxiliary probabilistic predictive model M . This approach was proposed recently by [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF]). The predictive model M defines a conditional probability distribution of target variables given the value of context variables: p(A|C) and it is induced (learned) from examples in E. Given M we say the case x is anomalous in target attributes A, if the probability p(A(x)|C(x)) for the model is small and falls below some threshold. In summary, stochastic relations in between the context and target attributes observed in examples E are incorporated into a probabilistic model M , which is in turn applied to example x. The anomaly is detected by evaluating the probability of target variable values for the x given the values of its context variables C(x) in model M .

To build a working anomaly detection algorithm, we need to provide methods for building a probabilistic model M from the dataset and methods for detecting the anomaly using the model.

Building a probabilistic model

Our conditional anomaly framework builds upon the existence of an underlying probabilistic model M that describes stochastic relations among context and target data attributes.

We consider two types of models to achieve this task: (1) parametric and (2) non-parametric models.

Parametric predictive models. In the parametric approach we assume a predictive model p(A|C) is defined using a small set of parameters Θ that reflect accurately the stochastic relation among the context and target attributes expressed in data E. Examples of parametric models are: a Bayesian belief network [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF], a Naive Bayes model [START_REF] Domingos | On the optimality of the simple bayesian classifier under zero-one loss[END_REF], Linear discriminant analysis [START_REF] Hastie | The Elements of Statistical Learning[END_REF]or a logistic regression model [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. In this work we focus on the Naive Bayes model that is used frequently in classification tasks. We adopt the Bayesian framework to learn the parameters of the model from data E and to support probabilistic inferences. In such a case the parameters M of the model are treated as random variables and are described in terms of a density function p(θ M ). The probability of an event is obtained by averaging over all possible parameter settings of the model M .

To incorporate the effect of examples E, p(θ M ) corresponds to the posterior p(θ M |E). The posterior is obtained via Bayes theorem:

p(θ M |E) = p(E|θ M )p(θ M )/p(E),
where p(θ M ) defines the prior for parameters θ M . To simplify the calculations we assume [START_REF] Heckerman | A tutorial on learning with bayesian networks[END_REF] (1) parameter independence and (2) conjugate priors. In such a case, the posterior follows the same distribution as the prior and updating reduces to updates of sufficient statistics. Similarly, many probabilistic calculations can be performed in the closed form.

Instance-specific models. In general, a predictive probabilistic model used for anomaly detection purposes can be of different complexity. However, if the dataset used to learn the model is relatively small, a more complex model may become hard to learn reliably. In such a case a simpler parametric model of P (A|C) with a smaller number of parameters may be preferred. Unfortunately, a simpler model may sacrifice some flexibility and its predictions may become biased towards the population of examples that occurs with a higher prior probability. To make accurate predictions for any instance we use instance-specific predictive methods and models [START_REF] Visweswaran | Instance-specific Bayesian model averaging for classification[END_REF][START_REF] Aha | Instancebased learning algorithms[END_REF].

Briefly, instance-based methods do not try to learn a universal predictive model for all possible instances, instead the model is optimized for every data instance x individually. To reflect this, we denote the predictive model for x as M x . The benefit of instance-based parametric models is that they can be fit more accurately to any data instance; the limitation is that the models must be trained only on the data that are relevant for x. Choosing the examples that are most relevant for training the instance-specific model is the bottleneck of the method. We discuss methods to achieve this later on.

Non-parametric predictive models. Non-parametric predictive models do not assume any compact parametrization of P (A|C). Instead, the model is defined directly on the dataset of examples E. A classic example of a non-parametric model is the k Nearest Neighbor (k-NN) classifier in which the predicted class of the instance is the majority vote of the classes of its k nearest neighbors.

Non-parametric models are instance specific by definition. For example, the k-NN classifier for instance x executes by finding k examples closest to x first and making the prediction afterwards. The problem of finding the k closest neighbors is the bottleneck of the method. Nonparametric models depend on the choice of examples closest to x, and the quality of these choices influences the quality of the model.

The anomaly detection approach applied in this work builds upon the model M x which defines the probability distribution P (A|C(x)) for x. But how to define a nonparametric predictive model M x ? The key here is to define the probability with which a neighbor example predicts the values of target attributes A for x. Intuitively, closer neighbors should contribute more and hence their prediction should come with a higher probability. To reflect this intuition [START_REF] Goldberger | Neighbourhood components analysis[END_REF] define the probability that a data example x ′ predicts x using the softmax model [START_REF] Mccullagh | Generalized Linear Models[END_REF]. In this model, the probability with which x ′ contributes to the prediction of x is proportional to:

exp(-||x -x ′ || 2 m )
where m is a distance metric reflecting the similarity of the samples. The softmax model normalizes this quantities so that their sum is 1.

The above definition of a non-parametric probabilistic predictive model expects a distance metric defining the similarity among examples. We return to the problem of distance metrics in the next section.

Anomaly detection

Multiple approaches can be used to make anomaly calls based on the probabilistic metric. Typically, they rely on a variety of thresholds. These include: absolute, relative or the k standard deviation thresholds. In our work, we build upon the absolute threshold test. In the absolute threshold test, the example x is anomalous if p(A(x)|C(x), M x ) falls below some fixed probability threshold p ε . Intuitively, if the probability of the target attributes A(x) for x is low with respect to the model M x and its other attributes C(x), then the value of the target attribute is anomalous. Note that the absolute threshold test relies only on the model M x and there is no direct comparison of predictive statistics for x and examples in E. However, if instance-based methods are used the most important examples in E are used to construct the model M x and hence their effect is reflected in the statistic. Exact match. Clearly, the best examples are the ones that exactly match the attributes C(x), of the target case x. How-ever, it is very likely that in real-world databases none or only few past cases match the target case exactly so there is no or very weak population support to draw any statistically sound anomaly conclusion.

Defining the similarity metric

Similarity-based match. One way to address the problem of insufficient population available through the exact match is to define a distance metric on the space of attributes C(x) that let us select examples closest to the target example x. The distance metric defines the proximity of any two cases in the dataset, and the k closest matches to the target case define the best population of size k. Different distance metrics are possible. An example is the generalized distance metric r 2 defined:

r 2 (x i , x j ) = (x i -x j ) T Γ -1 (x i -x j ),
(1) where Γ -1 is a matrix that weights attributes of patient cases proportionally to their importance. Different weights lead to a different distance metric. For example, if Γ is the identity matrix I, the equation defines the Euclidean distance of x i relative to x j . The Mahalanobis distance [START_REF] Mahalanobis | On the generalized distance in statistics[END_REF]) is obtained from 1 by choosing Γ to be the population covariance matrix Σ which lets us incorporate the dependencies among the attributes.

The Euclidean and Mahalanobis metrics are standard off-shelf distance metrics often applied in many learning tasks. However they come with many deficiencies. The Euclidean metric ignores feature correlates which leads to 'double-counting' when defining the distance in between the points. The Mahalanobis distance resolves this problem by reweighting the attributes according to their covariances. Nevertherless, the major deficiency of both Mahalanobis and Euclidean metrics is that they may not properly determine the relevance of an attribute for predicting the target attributes.

The relevance of context attributes for anomaly detection is determined by their influence on target attributes A. Intuitively, a context attribute is relevant for the predictive model if is able to predict changes in values of target attributes A. To incorporate the relevance aspect of the problem into the metric we adapt (learn) the parameters of the generalized distance metric with the help of examples in the dataset E.

Metric-learning

The problem of distance metric learning in context of classification tasks has been studied by [START_REF] Goldberger | Neighbourhood components analysis[END_REF]) and [START_REF] Bar-Hillel | Learning a mahalanobis metric from equivalence constraints[END_REF]. We adapt these metric learning methods to support probabilistic anomaly detection. In the following we briefly summarize the two methods. [START_REF] Goldberger | Neighbourhood components analysis[END_REF]) explores the learning of the metric in context of the nearest neighbor classification. They learn a generalized metric:

d 2 (x 1 , x 2 ) = (x 1 -x 2 ) T Q(x 1 -x 2 ) = (x 1 -x 2 ) T A T A(x 1 -x 2 ) = (Ax 1 -Ax 2 ) T (Ax 1 -Ax 2 )
by directly learning its corresponding linear transformation A. They introduce a new optimization criterion (NCA), that is, as argued by the authors, more suitable for the nearestneighbor classification purposes. The criterion is based on a new, probabilistic version of the cost function for the leaveone-out classification error in the k-NN framework. Each point i can now select any other point j with some probability p ij defined as softmax function over distances in the transformed space:

p ij = exp(-||Ax i -Ax j || 2 ) k =i exp((-||Ax k -Ax j || 2 )
A linear transformation A is then sought to maximize the expected number of correctly classified cases (with k-NN):

arg max A g(A) = arg max A i j∈Ci p ij
where C i is the set of cases that belong to the same class as i.

Intuitively, the criterion aims to learn a generalized distance metric by shrinking the distance between similar points to zero, and expanding the distance between dissimilar points to infinity. The algorithm and the metric it generates was shown to outperform other metrics for a number of learning problems. The method climbs the gradient of g(A), which is (x ij being x i -x j ): [START_REF] Bar-Hillel | Learning a mahalanobis metric from equivalence constraints[END_REF] and [START_REF] Shental | Adjustment learning and relevant component analysis[END_REF]) define a different optimization criterion based on the mutual information. The advantage of their method (relevant component analysis -RCA) is the existence of the closed form (efficient) solution. Briefly, under the mutual information criterion, the class information is incorporated and optimized by computing the averages of class covariance matrices. The resulting matrix is obtained by

∂g ∂A = 2A i   p i k p ik x ik x T ik - j∈Ci p ij x ij x T ij   (Bar
Σ RCA = k i=1 Σi A = Σ -1 2 (2)
where Σi sample covariance matrix of class i and A is the resulting transformation for the data. The disadvantage of the method is that it assumes Gaussian distribution for the classes.

Experimental evaluation

We study our metric-learning methods and compare them to alternative methods on the problem of identification of anomalous patient-management decisions for patients with community acquired pneumonia. The data used in the experiment come from the Pneumonia PORT dataset [START_REF] Kapoor | Assessment of the variantion and outcomes of pneumonia: Pneumonia patient outcomes research team (PORT) final report[END_REF][START_REF] Fine | A prediction rule to identify lowrisk patients with community-acquired pneumonia[END_REF]. The Pneumonia PORT dataset is based on the study conducted from October 1991 to March 1994 on 2287 patients with community-acquired pneumonia from three geographical locations at five medical institutions. The original PORT data were analyzed by [START_REF] Fine | A prediction rule to identify lowrisk patients with community-acquired pneumonia[END_REF], who derived a prediction rule with 30-day hospital mortality rate as the outcome. The authors developed a logistic regression model, which helped to identify 20 attributes that contribute the most to the mortality rate of pneumonia. To explore the anomaly detection methods, we have experimented with a simpler version of the PORT dataset that records, for every patient, only the attributes identified by Fine's study [START_REF] Fine | A prediction rule to identify lowrisk patients with community-acquired pneumonia[END_REF]. The attributes are summarized in Figure 1. All attributes are binary with true / false (positive / negative) values.

Our objective was to detect unusual admission decisions (treat the patient at home versus in the hospital) which are captured by the variable 'Hospitalization'.

Study design

To study the performance of our anomaly detection methods, we used 100 patient cases (out of a total of 2287 of cases). The cases picked for the study consisted of 21 cases that were found anomalous according to a simple Naive Bayes detector (with detection threshold 0.05) that was trained on all cases in the database. The remaining 79 cases were selected randomly from the rest of the database. Each of the 100 cases was then evaluated independently by a panel of three physicians. The physicians were asked whether they agree with the hospitalization decision or not. Using panel's answers, the admission decision was labeled as anomalous when (1) at least two physicians disagreed with the actual admission decision that was taken for a given patient case or (2) all three indicated they were unsure (gray area) about the appropriateness of the management decision. Out of 100 cases, the panel judged 23 as anomalous hospitalization d ecisions; 77 patient cases were labeled as not being anomalous. The assessment of 100 cases by the panel represented the correct assessment of unusual hospitalization decisions.

Experiments

All the experiments followed the leave-one-out scheme. That is, for each example in the dataset of 100 patient cases evaluated by the panel, we first learn the metric. Next, we identified the cases in E most similar to it with respect to that metric. The cases chosen were either the closest 40 cases, or all the other cases (2286) in the dataset. We then learned the NB model and calculated the posterior probability of the decision. Alternatively, we calculated the probability of the decision using the softmax model and the learned metric.

The target example was declared anomalous if its posterior probability value fell below the detection threshold. The anomaly calls made by our algorithms were compared to the assessment of the panel and resulting statistics (sensitivity, specificity) were calculated. To gain insight on the overall performance of each method we varied its detection threshold and calculated corresponding receiver operating characteristic (ROC). For the hospital deployment no all thresholds are acceptable. Consequently, for the evaluation we selected only that part of the ROC curve, that corresponds to specificity equal or greater than 95% (see Figure 2). The 95% specificity limit means that at most 1 in 20 normal cases analyzed may yield a false alarm.

Results and Discussion

Table 1 shows the ROC statistics for the feasible detection range. For the softmax model, the NCA metric outperformed all other methods, whether it was using all cases (patients) or just the closest 40. We ascribe it to the fact, that NCA uses class information to weigh the features. The only other method that used class information was RCA. However, RCA uses class information only to consider (and average) covariance matrices for each class separately. Therefore, it still treats all features within the class the same way as the Mahalanobis metric, assuming the same relevance of all features. Comparing the global (all other patients) and local (closest 40 patients), local did always better: Close patients let us fit better the predictive model to the target patient, while taking all samples into the consideration biases the population. The local methods were also better for the Naive Bayes model. They were also more robust with respect to the metric. The intuition behind this result is that when using NB model, all cases are treated the same way, the metric was only used to select them. On the other hand, in softmax model, the distance from the case in hand does matter and the method treats closer patients with a higher weight. Accordingly, it is more sensitive to the metric changes.

Figure 2 shows the ROC curve for the best method in corresponds to the performance of 6.66 correct alarms in 10 alarms for 100 evaluated patients. However, we note that the prior for the evaluation dataset was biased towards anomalies. A rough correction using only anomalies that were randomly selected from the full database yields approximately 1 correct in 4 alarms, which is still very encouraging performance.

Conclusions

Summing up, our conditional anomaly detection is a very promising methodology for detecting unusual events such as network attacks or medical errors. We have demonstrated its potential by exploring and analyzing patient-management decisions for a dataset of patients suffering from pneumonia. The advantage of the anomaly detection approach over knowledge-based error detection approaches is that the method is evidence-based, and hence requires no or minimum input from the domain expert.

Despite initial encouraging results, our current approach can be further refined and extended. For example, instancebased (local) models tested in this paper always used a fixed number of 40 closest patients (or more, if the distances were the same). However, the patient's neighborhood and its size depend on the patient and data available in the database. We plan to address the problem by developing methods that are able to automatically identify and select only patients that are close enough for the case in hand.

  Parametric instance-based models are sensitive of examples used to train them. Similarly, the instance-based nonparametric models are sensitive to examples incorporated into the model. The key question is what examples from E should be used for training or defining the instance-specific predictive model M x .

Figure 1 :

 1 Figure 1: Attributes from the Pneumonia PORT dataset used in the anomaly detection study.

Figure 2 :

 2 Figure 2: An example of the ROC curve for the method that performed the best on the pneumonia dataset. The statistic of interest is the leftmost region of the ROC curve and its area.

Table 1 .

 1 The area of interest is bounded by the values [0, 0.13], [0.03, 0.30], and [0.05, 0.33]. The point [0.05 0.33]

	model selection
	non-parametric	global #cases	area
	NCA	softmax	2286 18.0 %
	Mahal	softmax	2286 12.2 %
	RCA	softmax	2286 11.6 %
	Euclidean	softmax	2286	8.0 %
	non-parametric	local #cases	area
	NCA	softmax	40 20.2 %
	Mahalanobis	softmax	40 15.0 %
	RCA	softmax	40 12.8 %
	Euclidean	softmax	40	8.0 %
	parametric	global #cases	area
	any	NB	2286 11.6 %
	parametric	local #cases	area
	NCA	NB	40 16.8 %
	Mahalanobis	NB	40 17.6 %
	RCA	NB	40 17.6 %
	Euclidean	NB	40 16.4 %

Table 1 :

 1 Area under the ROC curve in the feasible range of 95% -100% specificity. Please note that the baseline value for the random choice is 2.5%.

Acknowledgements

The research presented in this paper was funded by the grant R21-LM009102-01A1 from the National Library of Medicine.


Mathias Rousset
email: mathias.rousset@inria.fr

Giovanni Samaey
email: giovanni.samaey@cs.kuleuven.be

Simulating individual-based models of bacterial chemotaxis with asymptotic variance reduction

Keywords: 35Q80, 92B05, 65C35 asymptotic variance reduction, bacterial chemotaxis, velocity-jump process

We discuss variance reduced simulations for an individual-based model of chemotaxis of bacteria with internal dynamics. The variance reduction is achieved via a coupling of this model with a simpler process in which the internal dynamics has been replaced by a direct gradient sensing of the chemoattractants concentrations. In the companion paper 26 , we have rigorously shown, using a pathwise probabilistic technique, that both processes converge towards the same advection-diffusion process in the diffusive asymptotics. In this work, a direct coupling is achieved between paths of individual bacteria simulated by both models, by using the same sets of random numbers in both simulations. This coupling is used to construct a hybrid scheme with reduced variance. We first compute a deterministic solution of the kinetic density description of the direct gradient sensing model; the deviations due to the presence of internal dynamics are then evaluated via the coupled individual-based simulations. We show that the resulting variance reduction is asymptotic, in the sense that, in the diffusive asymptotics, the difference between the two processes has a variance which vanishes according to the small parameter.

Introduction

The motion of flagellated bacteria can be modeled as a sequence of run phases, during which a bacterium moves in a straight line at constant speed. Between two run phases, the bacterium changes direction in a tumble phase, which takes much less time than the run phase and acts as a reorientation. To bias movement towards regions with high concentration of chemoattractant, the bacterium adjusts its turning rate by increasing, resp. decreasing, the probability of tumbling when moving in an unfavorable, resp. favorable, direction [START_REF] Alt | Orientation of cells migrating in a chemotactic gradient[END_REF][START_REF] Stock | A nonlinear stimulus-response relation in bacterial chemotaxis[END_REF] . Since many species are unable to sense chemoattractant gradients reliably due to their small size, this adjustment is often done via an intracellular mechanism that allows the bacterium to retain information on the history of the chemoattractant concentrations along its path [START_REF] Bren | How signals are heard during bacterial chemotaxis: Protein-protein interactions in sensory signal propagation[END_REF] . The resulting model, which will be called the "internal state" or "fine-scale" model in this text, can be formulated as a velocity-jump process, combined with an ordinary differential equation (ODE) that describes the evolution of an internal state that incorporates this memory effect [START_REF] Erban | From signal transduction to spatial pattern formation in e-coli: A paradigm for multiscale modeling in biology[END_REF] .

The probability density distribution of the velocity-jump process evolves according to a kinetic equation, in which the internal variables appear as additional dimensions. A direct deterministic simulation of this equation is therefore prohibitively expensive, and one needs to resort to a stochastic particle method. Unfortunately, a direct fine-scale simulation using stochastic particles presents a large statistical variance, even in the diffusive asymptotic regime, where the behavior of the bacterial density is known explicitly to satisfy a Keller-Segel advection-diffusion equation. Consequently, it is difficult to assess accurately how the solutions of the fine-scale model differ from their advection-diffusion limit in intermediate regimes. We refer to [START_REF] Horstman | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences I[END_REF][START_REF] Horstman | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences II[END_REF][START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF] for numerous historical references on the Keller-Segel equation, to [START_REF] Erban | From individual to collective behavior in bacterial chemotaxis[END_REF][START_REF] Erban | From signal transduction to spatial pattern formation in e-coli: A paradigm for multiscale modeling in biology[END_REF] for formal derivations (based on moment closures) of convergence of the the velocityjump process with internal state to the Keller-Segel equation, and to the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] for a proof of this convergence based on probabilistic arguments.

In this paper, we propose and analyze a numerical method to simulate individual-based models for chemotaxis of bacteria with internal dynamics with reduced variance. The variance reduction is based on a coupling technique (control variate): the main idea is to simultaneously simulate, using the same random numbers, a simpler, "coarse" process where the internal dynamics is replaced by a direct "gradient sensing" mechanism (see [START_REF] Alt | Biased random-walk models for chemotaxis and related diffusion approximations[END_REF][START_REF] Hg Othmer | Models of dispersal in biological-systems[END_REF][START_REF] Patlak | Random walk with persistence and external bias[END_REF] for references on such gradient sensing models). The probability density of the latter satisfies a kinetic equation without the additional dimensions of the internal state, and converges to a similar advection-diffusion limit, see e.g. [START_REF] Facc Chalub | Kinetic models for chemotaxis and their drift-diffusion limits[END_REF][START_REF] Hillen | The diffusion limit of transport equations derived from velocity-jump processes[END_REF][START_REF] Hg Othmer | The diffusion limit of transport equations ii: Chemotaxis equations[END_REF][START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] .

More precisely, we consider the coupling of two velocity-jump processes with exactly the same advection-diffusion limit (for both processes, a proof of the latter convergence has been obtained in the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] ). The first one is the process with internal dynamics described by the system of equations (2.4), with jump rates satisfying (2.3) and internal dynamics satisfying the two main assumptions (2.7) and (2.9) as well as other technical assumptions contained in Section 2.2. The second one is the process with direct gradient sensing (2.12), with jump rates satisfying (2.14) with (2.15). The same random numbers are used for the two processes in the definition of both the jump times and the velocity directions. The main contributions of the present paper are then twofold :

• From a numerical point of view, we couple two systems of many particles consisting of different realizations of the fine-scale (internal state) and coarse (gradient sensing) processes, simulated with similar discretization schemes (see Section 3.1) and the same random numbers. Additionally, for the coarse model, also the continuum description is simulated on a spatial grid. Evolution of the fine-scale model is then evaluated on the grid at hand by adding the difference between the two coupled particle descriptions to the evolution of the coarse continuum description. This can be seen as using the coarse process as a control variate. As coupling is lost over time, a regular reinitialization of the two coupled particle systems is included in Model Numerics Internal state model-Individual description Monte-Carlo method Gradient sensing model-Individual description Monte-Carlo method (coupled) Gradient sensing model-Density description Grid method

Table 1. Models and numerical methods used in the variance reduction algorithm.

the algorithm. The asymptotic variance reduction method can be depicted using Table 1 (see also Figure 1 in Section 3). • From an analysis point of view, we rigorously prove, using probabilistic arguments, L p -bounds on the position difference at diffusive times of both processes in terms of appropriate powers of the scaling parameter (Theorem 4.7). This analysis shows that the numerical variance reduction obtained by the method is asymptotic (the remaining statistical variance vanishes with the small scaling parameter), and provides upper bounds on the rate of the latter convergence. We argue in Section 4.3 that this bound should be sharp in relevant regimes.

The idea of asymptotic variance reduction is a general and very recent idea in scientific computing that appears when using hybrid Monte-Carlo/PDE (partial differential equation) methods. While the use of control variates is fairly common in Monte Carlo simulation, the only explicit attempt to develop asymptotic variance reduction in hybrid methods, up to our knowledge, can be found in [START_REF] Dimarco | Hybrid multiscale methods II. Kinetic equations[END_REF] and related papers in the context of the Boltzmann equation, and is based on an importance sampling approach (see also [START_REF] Baker | Variance reduction for monte carlo solutions of the boltzmann equation[END_REF] and related papers).

This paper is organized as follows. In Section 2, we briefly review the mathematical models that we will consider, and summarize the results on their advectiondiffusion limit that were obtained in the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] . In Section 3, we proceed to describe the coupled numerical method in detail. Section 4 is dedicated to a detailed analysis of the variance reduction of the coupled method. The main result is the following: up to a time shift perturbation, the variance of the difference between the two coupled process scales according to the small parameter of the diffusive asymptotics. We proceed with numerical illustrations in Sections 5 and 6, and conclude in Section 7 with an outlook to future research.

Particle-based models for bacterial chemotaxis

In this paper, we directly present the models of interest in a nondimensional form, incorporating the appropriate space and time scales. For a discussion on the chosen scaling, we refer to the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] . Variance reduced simulation of chemotaxis

Model with internal state

We consider bacteria that are sensitive to the concentration of m chemoattractants (ρ i (x)) m i=1 , with ρ i (x) ≥ 0 for x ∈ R d . While we do not consider time dependence of chemoattractant via production or consumption by the bacteria, a generalization to this situation is straightforward, at least for the definition of the models and the numerical method. Bacteria move with a constant speed v (run), and change direction at random instances in time (tumble), in an attempt to move towards regions with high chemoattractant concentrations. As in [START_REF] Erban | From individual to collective behavior in bacterial chemotaxis[END_REF] , we describe this behavior by a velocity-jump process driven by some internal state y ∈ Y ⊂ R n of each individual bacterium. The internal state models the memory of the bacterium and is subject to an evolution mechanism attracted by a function ψ : R m → R n of the chemoattractants concentrations,

S(x) := ψ(ρ 1 (x),...,ρ m (x)) ∈ R n ,
where x is the present position of the bacterium. S is assumed to be smooth with bounded derivatives up to order 2. We refer to Section 2.2 in 26 for some typical choices; in the concrete example in this paper, we choose m = n = 1, and S(x) = ρ 1 (x). By convention, the gradient of S(x) : R d → R n is a matrix with dimension

∇S(x) ∈ R n×d .
We then denote the evolution of each individual bacterium position by t → X t , with normalized velocity

dX t dt = ǫV t , V t ∈ V = S d-1 ,
with S d-1 the unit sphere in R d . Hence, V t represents the direction and the parameter ǫ represents the size of the velocity. The evolution of the internal state is denoted by t → Y t . The internal state adapts to the local chemoattractant concentration through an ordinary differential equation (ODE),

dY t dt = F ǫ (Y t ,S(X t )), (2.1)
which is required to have a unique fixed point y * = S(x * ) for every fixed value x * ∈ R d . We also introduce the deviations from equilibrium z = S(x) -y. The evolution of these deviations is denoted as

t → Z t = S(X t ) -Y t .
The velocity of each bacterium is switched at random jump times (T n ) n≥1 that are generated via a Poisson process with a time dependent rate given by λ(Z t ), where z → λ(z) is a smooth function satisfying

0 < λ min ≤ λ(z) ≤ λ max , (2.2)
as well as,

λ(z) = λ 0 -b T z + c λ O |z| k , (2.3)
with b ∈ R n , k ≥ 2, and c λ > 0 is used to keep track of the non-linearity in the analysis. The new velocity at time T n is generated at random according to a centered probability distribution M(dv) with v M(dv) = 0, typically

M(dv) = σ S d-1 (dv),
where σ S d-1 is the uniform distribution on the unit sphere.

The resulting fine-scale stochastic evolution of a bacterium is then described by a left continuous with right limits (lcrl) process t → (X t ,V t ,Y t ), which satisfies the following differential velocity-jump equation :

                       dX t dt = ǫV t , dY t dt = F ǫ (Y t ,S(X t )), Tn+1 Tn λ(Z t )dt = θ n+1 , with Z t := S(X t ) -Y t , V t = V n for t ∈ [T n ,T n+1 ] , (2.4)
with initial condition X 0 ,V 0 ∈ R d , Y 0 ∈ R n and T 0 = 0. In (2.4), (θ n ) n≥1 denote i.i.d. random variables with normalized exponential distribution, and (V n ) n≥1 denote i.i.d. random variables with distribution M(dv).

Example 2.1. For concreteness, we provide a specific example, adapted from 11 , which will also be used later on to illustrate our results numerically. We consider m = 1, i.e., there is only one chemoattractant S(x) and the internal dynamics (2.1) reduces to a scalar equation

dy(t) dt = S(x) -y τ = z τ . (2.5)
For the turning rate z → λ(z), we choose the following nonlinear strictly decreasing smooth function

λ(z) = 2λ 0 1 2 - 1 π arctan π 2λ 0 z . (2.6)

Asymptotic regimes and linearizations

In the adimensional models (2.4)-(2.12), a small parameter ǫ is introduced, which can be interpreted as the ratio of the typical time of random change of velocity direction to the typical time associated with chemoattractant variations as seen by the bacteria. In the present section, we present the main scaling assumptions on the internal dynamics that are necessary to consider the asymptotic regime ǫ ≪ 1, as well as the corresponding notation. Most of these assumptions were introduced in the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] , to which we refer for more details and motivation. Variance reduced simulation of chemotaxis Remark 2.1 (Notation). Throughout the text, the Landau symbol O denotes a deterministic and globally Lipschitz function satisfying O(0) = 0. Its precise value may vary from line to line, may depend on all parameters of the model, but its Lipschitz constant is uniform in ǫ. In the same way, we will denote generically by C > 0 a deterministic constant that may depend on all the parameters of the model except ǫ.

We assume that the ODE (2.1) driving the internal state is well approximated by a near equilibrium evolution equation in following sense : Assumption 2.2. We have

F ǫ (y,s) = -τ ǫ -1 (y -s) + ǫ 1-δ c F O(|s -y| 2 ), (2.7)
where δ > 0 and τ ǫ ∈ R n×n is an invertible constant matrix.

The constant c F is used to keep track of the dependance on the non-linearity of F ǫ in the analysis. For technical reasons that are specific to the analysis of the time shift in the coupling (Lemma 4.6), we need the following assumption on the scale of the non-linearity (this assumption was not necessary in the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] ): Assumption 2.3. Assume that |y -s| = O(ǫ δ ), then there is a γ > 0 such that:

|(y -s) -τ ǫ F ǫ (y,s)| = c F O(ǫ 1+γ ).
(2.8)

We now formulate the necessary assumptions on the scale of the matrix τ ǫ : Assumption 2.4. There is a constant C > 0 such that for any t ≥ 0, one has e -tτǫ -1 ≤ Ce -tǫ 1-δ /C . (2.9) Assumption 2.5. There is a constant C > 0 such that

sup t≥0 tτ ǫ -1 e -tτǫ -1 ≤ Cǫ -1 .
Assumption 2.4 is used to ensure exponential convergence of the linear ODE with time scale at least of order ǫ 1-δ . Assumption 2.5 is necessary for technical reasons in the proof of Lemma 4.5 (which is given in [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] ). Note that if τ ǫ -1 is symmetric and strictly positive with a lower bound of order O(ǫ 1-δ ) on the spectrum, then Assumption 2.4 and Assumption 2.5 are satisfied. This implies that the slowest time scale of the ODE (2.1) is at worse of order ǫ 1-δ . Indeed, integration over time in Assumption 2.4 implies τ ǫ ≤ Cǫ δ-1 . Moreover, we assume that the solution of the ODE driving the internal state in (2.4) τ ǫ -1 e -τǫ -1 t Z t = O(1).

τ ǫ -1 e -τǫ -1 t (S t -Y t ) = O (1), (2.10)
Finally, we will make use of the following notation. Denoting c S = ∇ 2 S ∞ , with ∇ 2 the Hessian, the error terms due to non-linearities of the problem will be handled through

nl(ǫ) := c F ǫ δ + c S ǫ + c λ ǫ kδ-1 = o(1), nl 2 (ǫ) := c F (ǫ δ + ǫ γ ) + c S (ǫ + ǫ δ ) + c λ ǫ kδ-1 = o(1).
(2.11)

Note that nl(ǫ) ≤ nl 2 (ǫ) and nl(ǫ) ∼ nl 2 (ǫ) if δ ≤ 1 and δ ≤ γ. Recall that k ≥ 2 is defined in (2.3).
Note that all the assumptions of the present section hold in the following case: (i) the ODE is linear (F ǫ is linear), (ii) the associated matrix τ ǫ is symmetric positive and satisfies

τ -1 ǫ ≥ Cǫ 1-δ ,
for some δ > 0, (iii) the initial internal state is close to equilibrium in the sense that τ ǫ -1 Z 0 = O(1). The assumptions of the present section can be seen as technical generalizations to non-linear ODEs with non-symmetric linear part. Variance reduced simulation of chemotaxis 2.3. Model with direct gradient sensing (control process)

We now turn to a simplified, "coarse" model, in which the internal process (2.1), and the corresponding state variables, are eliminated. Instead, the turning rate depends directly on the chemoattractant gradient. This process will be called the control process since it will be used in Section 3 as a control variate to perform variance reduced simulations of (2.4).

The control process is a Markov process in position-velocity variables t → (X c t ,V c t ), which evolves according to the following differential velocity-jump equations :

             dX c t dt = ǫV c t , T c n+1 T c n λ c ǫ (X c t ,V c t )dt = θ n+1 , V c t = V n for t ∈ [T c n ,T c n+1 ] , (2.12)
with initial condition X 0 ,V 0 ∈ R d . In (2.12), (θ n ) n≥1 denote i.i.d. random variables with normalized exponential distribution, and (V n ) n≥1 denote i.i.d. random variables with distribution M(dv). The turning rate of the control process is assumed to satisfy

0 < λ min ≤ λ c ǫ (x,v) ≤ λ max , (2.13)
as well as

λ c ǫ (x,v) := λ 0 -ǫ A T ǫ (x)v + O(ǫ 2 ). (2.14)
Typically, λ c ǫ (x,v) is a function of ∇S(x), so that the model (2.14) may describe a large bacterium that is able to directly sense chemoattractant gradients. When m = n = 1, and the turning rate (2.14) is proportional to ∇S(x)v ∈ R, it can be interpreted as follows: the rate at which a bacterium will change its velocity direction depends on the alignment of the velocity with the gradient of the chemoattractant concentration ∇S(x), resulting in a transport towards areas with higher chemoattractant concentrations. In the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] , it is shown that, when the vector field A ǫ (x) ∈ R d is given by the formula

A ǫ (x) = b T τ ǫ Id+λ 0 τ ǫ ∇S(x), (2.15)
asymptotic consistency with the internal state model is obtained in the limit when ǫ → 0 (see Section 2.4).

Kinetic formulation and advection-diffusion limits

We now turn to the kinetic description of the probability distribution of the processes introduced above. The probability distribution density of the fine-scale process with internal state at time t with respect to the measure dxM(dv)dy is denoted as p(x,v,y,t), suppressing the dependence on ǫ for notational convenience, and evolves according to the Kolomogorov forward evolution equation (or master equation). In the present context, the latter is the following kinetic equation

∂ t p + ǫv • ∇ x p + div y (F ǫ (x,y)p) = λ(S(x) -y)(R(p) -p), (2.16)
where

R(p) := v∈S d-1 p(•,v,•)M(dv)
is the operator integrating velocities with respect to M, and F (x,y) and λ(z) are defined as in Section 2.1. Similarly, the distribution density of the control process is denoted as p c (x,v,t), and evolves according to the kinetic/master equation:

∂ t p c + ǫv • ∇ x p c = (R(λ c ǫ p) -λ c ǫ p), (2.17)
with λ c ǫ defined as in (2.14). The reader is referred to [START_REF] Ethier | Markov Processes: Characterization and Convergence[END_REF] for the derivation of master equations associated to Markov jump processes.

In the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] , we show, using probabilistic arguments, that, in the limit ǫ → 0, both the equation for the control process (2.17) and the equation for the process with internal state (2.16) converge to an advection-diffusion limit on diffusive time scales. Convergence is to be understood pathwise, in the sense of convergence of probability distribution on paths endowed with the uniform topology.

Denote diffusive times by t := tǫ 2 , and the processes considered on diffusive time scales as

X ǫ t := Xt /ǫ 2 , X c,ǫ t := X c t/ǫ 2 .
For the control process, we then have the following theorem :

Theorem 2.1. For ǫ → 0, the process t → X c,ǫ t , solution of (2.12), converges towards an advection-diffusion process, satisfying the stochastic differential equation (SDE)

dX c,0 t = DA 0 (X c,0 t ) λ 0 d t + 2D λ 0 1/2 dWt , (2.18)
where t → Wt is a standard Brownian motion, the parameters λ 0 and A 0 := lim ǫ→0 A ǫ originate from the turning rate (2.14), and the diffusion matrix is given by the covariance of the Maxwellian distribution :

D = S d-1 v ⊗ v M(dv) ∈ R d×d . (2.19)
In particular, this result implies that, at the level of the Kolomogorov/master evolution equation, the evolution of the position bacterial density,

n c,ǫ (x, t) := n c (x, t/ǫ 2 ) := V p c (x,v, t/ǫ 2 )M(dv)
(2.20) Variance reduced simulation of chemotaxis converges to the advection-diffusion equation

∂tn c,0 = 1 λ 0 div x D∇ x n c,0 -DA 0 (x)n c,0 (2.21)
on diffusive time scales as ǫ → 0.

In the same way, a standard probabilistic diffusion approximation argument can be used to derive the pathwise diffusive limit of the process with internal state (2.4) : Theorem 2.2. For ǫ → 0, the process t → X ǫ t , solution of (2.4), converges towards an advection-diffusion process, satisfying the stochastic differential equation (SDE) (2.18), where A 0 originates from

A 0 (x) = b T lim ǫ→0 τ ǫ λ 0 τ ǫ + Id ∇S(x), (2.22)
in which, b, τ ǫ , and λ 0 were introduced in (2.7)-( 2.3) as parameters of the process with internal state, and Id ∈ R n×n is the identity matrix. Again, the diffusion matrix D is given by the covariance of the Maxwellian distribution (2.19).

Also here, convergence needs to be understood in terms of probability distribution on paths endowed with the uniform topology.

Introducing the bacterial density of the process with internal state as

n(x,t) = Y V p(x,v,y,t)M(dv)dy, (2.23)
this implies that the evolution of n converges to (2.21) on diffusive time scales in the limit of ǫ → 0.

Numerical method

To simulate the process with internal state, solving the kinetic equation (2.16) over diffusive time scales can be cumbersome, due to the additional dimensions associated with the internal state. The alternative is to use to stochastic particles. However, a particle-based simulation of equation (2.16) is subject to a large statistical variance of the order O(N -1/2 ), where N is the number of simulated particles.

The asymptotic analysis shows that the position bacterial density approaches an advection-diffusion limit (2.21) when ǫ → 0. Consequently, to accurately assess the deviations of the process with internal state (2.4) as compared to its advectiondiffusion limit (for small but non-vanishing (intermediate) values of ǫ), the required number of particles needs to increase substantially with decreasing ǫ, which may become prohibitive from a computational point of view.

In this section, we therefore propose a hybrid method, based on the principle of control variates, that couples the process with internal dynamics with the control process, which is simulated simultaneously using a grid-based method. We first describe the variance reduction technique (Section 3.1). The analysis in Section 4 will reveal that the variance reduction is asymptotic, in the sense that the variance vanishes in the diffusion limit. To ensure this asymptotic variance reduction during actual simulations, one needs to ensure that the time discretization preserves the diffusion limits of the time-continuous process. An appropriate time discretization is discussed in Section 3.2.

Coupling and asymptotic variance reduction

The proposed variance reduction technique is based on the introduction of a control variate that exploits a coupling between the process with internal state and the control process. While the idea of control variates for Monte Carlo simulation is already well known, see e.g. [START_REF] Kloeden | Numerical solution of stochastic differential equations[END_REF] and references therein, the coupling that is proposed here is particular (we call it asymptotic), since the difference of the coupled processes, and hence the variance, vanishes with an estimable rate in the diffusion limit ǫ → 0, as will be shown in Section 4.

The control variates

Let us first assume that we are able to compute the exact solution of the kinetic equation for the control process, (2.17), with infinite precision in space and time.

The algorithm of asymptotic variance reduction is based on a coupling between an ensemble of realizations evolving according to the process with internal state (2.4), denoted as

X i t ,V i t ,Y i t N i=1
, and an ensemble of realization of the control process (2.12), denoted as

X i,c t ,V i,c t N i=1
.

We denote the empirical measure of the particles with internal state in positionvelocity space as

µ N t (x,v) = 1 N N i=1 δ X i t/ǫ 2 ,V i t/ǫ 2 ,
and, correspondingly, the empirical measure of the control particles as

µ c,N t = 1 N N i=1 δ X c,i t/ǫ 2 ,V c,i t/ǫ 2 .
A coupling between the two ensembles is obtained by ensuring that both simulations use the same random numbers (θ n ) n≥1 and (V n ) n≥0 , which results in a strong correlation between (X i t ,V i t ) and (X i,c t ,V c,i t ) for each realization. Simultaneously, the kinetic equation for the control process (2.17) is also solved using a deterministic method (which, for now, is assumed to be exact). We formally denote the corresponding semi-group evolution as

e tL c , with L c (p c ) = -ǫv • ∇ x p c + (R(λ c ǫ p c ) -λ c ǫ p c ).
Variance reduced simulation of chemotaxis

Besides the two particle measures µ N t and µ c,N t , we denote by µ N t the variance reduced measure, which will be defined by the algorithm below. Since, with increasing diffusive time, the variance of the algorithm increases due to a loss of coupling between the particles with internal state and the control particles, the variance reduced algorithm will also make use of a reinitialization time step δt ri , which is defined on the diffusive time scale. The corresponding time instances are denoted as tn = nδt ri on the diffusive time scale, or equivalently, on the original time scale as

t n = nδt ri /ǫ 2 .
Starting from an initial probability measure µ 0 at time t = 0, we sample µ 0 to obtain the ensemble

X i t ,V i t ,Y i t N i=1
, corresponding to µ N 0 , and then set µ c,N 0 := µ N 0 , i.e., X i,c 0 = X i 0 and V i,c 0 = V i 0 for all i = 1,...,N . Furthermore, we set the variance reduced estimator as µ N 0 := µ 0 = E(µ N 0 ). We then use the following algorithm to advance from tn to tn+1 , (see also Figure ( 1)) : Algorithm 3.1. At time t n , we have that the particle measure µ c,N tn = µ N tn , and the variance reduced measure is given by µ N tn . To advance from time tn to tn+1 , we perform the following steps :

• Evolve the particles X i t ,V i t ,Y i t N i=1 from t n to t n+1 , according to (2.4), • Evolve the particles X i,c t ,V i,c t N i=1
according to (2.12), using the same random numbers as for the process with internal state, • Compute the variance reduced evolution

µ N tn+1 = µ N tn e δtri/ǫ 2 L c + µ N tn+1 -µ c,N t- n+1 . (3.1)
• Reinitialize the control particles by setting

X i,c tn+1 = X i tn+1 , V i,c tn+1 = V i tn+1 , i = 1,...,N,
i.e., we set the state of the control particles to be identical to the state of the particles with internal state.

In formula (3.1), we use the symbol t-n+1 to emphasize that the involved particle positions and velocities are those obtained before the reinitialization. An easy computation shows that the algorithm is unbiased in the sense that for any n ≥ 0,

E µ N tn = E µ N tn ,
since the particles with internal dynamics are unaffected by the reinitialization, and, additionally,

Eµ c,N tn+1 = Eµ N tn e δtri/ǫ 2 L c .
Moreover, the variance is controlled by the coupling between the two processes. Indeed, using the independence of the random numbers between two steps of Algorithm 3.1, and introducing ϕ as a position and velocity dependent test function), The dotted line represent the coupled evolution of N bacteria with gradient sensing, subject to regular reinitializations. The dashed-dotted line is computed according to a deterministic method simulating the density of the model with gradient sensing, and subject to regular reinitializations. The solid line is the variance reduced simulation of the internal state dynamics, and is computed by adding the difference between the particle computation with internal state, and the particle simulation with gradient sensing to the deterministic gradient sensing simulation. At each reinitialization step, the two simulations (deterministic and particles) of the gradient sensing dynamics are reinitialized to the values of their internal state simulation counterpart (as represented by the arrows).

we get (according to the main result of Theorem 4.7)

var(µ N tn (ϕ)) = n k=1 E µ N (ϕ) -µ c,N t- k (ϕ) 2 ≤ n k=1 ∇ϕ ∞ N E Xt k /ǫ 2 -X c tk /ǫ 2 2
, Variance reduced simulation of chemotaxis and thus

var(µ N tn (ϕ)) ≤ Cn ǫ + ǫ δ + nl 2 (ǫ) N , (3.2)
where in the last line, C is independant of n, ǫ, and N .

In some generic situations (see Section 4.3), we can argue that the statistical error in Algorithm 3.1 coming from the coupling is "sharp" with respect to the order in ǫ. This means that the difference between the probability distribution of the model with internal state and the probability distribution of the model with gradient sensing is of the same order. This would imply that, with the asymptotic variance reduction technique, one is able to reliably assess the true deviation of the process with internal variables from the control process using a number of particles N that is independent of ǫ.

Time-space discretization of the kinetic equation

For both processes, the probability density distribution of the position and velocity can then be computed via binning in a histogram, or via standard kernel density estimation 27,28 using a kernel K h , in which the chosen bandwidth h can be based on the grid used for the deterministic simulation and on the data, for instance using the Silverman heuristic [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF] . This yields

pN (x,v,t) = 1 N h N i=1 K h (x -X i t ,v -V i t ), (3.3)
and similarly

pc N (x,v,t) = 1 N h N i=1 K h (x -X c,i t ,v -V c,i t ). (3.4)
The important point is that the solution p c (x,v,t) for the control process may be accurately approximated with a deterministic grid-based method, which ensures that

p c (x,v, t/ǫ 2 ) -p c grid (x,v, t/ǫ 2 ) = O(δx l ) + O(δt k ),
for some integers k,l ≥ 1.

Remark 3.1. In this text, we will perform the simulations in one space dimension using a third-order upwind-biased scheme, and perform time integration using the standard fourth order Runge-Kutta method, i.e., l = 3 and k = 4.

Then, the unbiased nature of the variance-reduced estimator is conserved up to O(δx l ) + O(δt k ) + O(h) discretization errors.

Time discretization of velocity-jump processes

When simulating equation (2.4), a time discretization error originates from the fact that the equation for the evolution of the internal state, and hence the evolution of the fluctuations Z t , is discretized in time. This results in an approximation of the jump times (T n ) n≥1 , and hence of X t . To retain the diffusion limits of the time-continuous process, special care is needed. We now briefly recall the time discretization procedure that was proposed and analyzed in the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] . For ease of exposition, we consider the scalar equation (2.5) for the internal state; generalization to nonlinear systems of equations is briefly discussed in [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] .

Consider first the linear turning rate (2.3). In that case, we define a numerical solution (X δt t ,V δt t ,Y δt t ) as follows. Between jumps, we discretize the simulation in steps of size δt and denote by (X δt n,k ,Z δt n,k ) the solution at t n,k = T δt n + kδt. The numerical solution for t ∈ [t n,k ,t n,k+1 ] is given by

X δt t = X δt n,k + ǫV n (t -t n,k ) Z δt t = exp(-(t -t n,k )τ ǫ -1 )Z δt n,k + ǫτ ǫ Id-exp(-(t -t n,k )τ ǫ -1 ) ∇S(X δt n,k )V n . (3.5)
We denote by K ≥ 0 the integer such that the simulated jump time T δt n+1 ∈ [t n,K ,t n,K+1 ]. To find T δt n+1 , we first approximate the integral

T δt n+1 T δt n λ(Z δt t )dt us- ing T δt n+1 T δt n λ(Z δt t )dt = K-1 k=0 t n,k+1 t n,k λ(Z δt t )dt + T δt n+1 t n,K λ(Z δt t )dt, (3.6)
and then compute :

t n,k+1 t n,k λ(Z δt t )dt = λ 0 δt -b T Id-e -δtτǫ -1 τ ǫ Z δt Tn -ǫb T δtτ ǫ -(Id-e -δtτǫ -1 )τ ǫ 2 ∇S(X δt n,k )V n . (3.7)
The jump time T δt n+1 can then be computed as the solution of

T δt n+1 t n,K λ(Z δt t )dt = θ n+1 - K-1 k=0 t n,k+1 t n,k λ(Z δt t )dt, (3.8)
using a Newton procedure. It is shown in [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] that the results on the diffusion limit (as outlined in Section 2.4) are not affected by the discretization. We now consider a general nonlinear turning rate. We again discretize in time to obtain the time-discrete solution (3.5). The jump time T δt n+1 is now computed by linearizing (2.6) in each time step,

λ δt Z δt t ,Z δt n,k = λ(Z δt n,k ) + dλ(Z δt n,k ) dz Z δt t -Z δt n,k .
(3.9) Variance reduced simulation of chemotaxis

One then replaces equation (3.6) by

T δt n+1 T δt n λ δt (Z δt t )dt = K-1 k=0 t n,k+1 t n,k λ δt (Z δt t ,Z δt n,k )dt + T δt n+1 t n,K λ δt (Z δt t ,Z δt n,K )dt, (3.10)
and proceeds in the same way as for the linear case. Also in this case, the diffusive limit is recovered in an exact fashion for the time-discretized process.

Remark 3.2 (Discretization of the control process). For the control process, there is no internal state. Hence, the only time dependence of turning rate λ c is due to the spatial variation of ∇S(x), which can be treated by discretizing the integral of λ in the same way as above.

Asymptotic variance reduction of the coupling

In this section, we show how the difference between the two coupled processes on diffusive time scales t = t/ǫ 2 behaves in the limit of ǫ → 0. We first recall some notation from the companion paper 26 (Section 4.1), after which we state and prove the main theorem (Section 4.2).

Notations and asymptotic estimates

The main theorem in Section 4.2 relies on following auxiliary definitions and lemmas that were given in the companion paper 26 : Definition 4.1. We denote by m : R → R n×n the function

m(t) := tτ ǫ -Id-e -tτǫ -1 τ ǫ 2 , (4.1)
whose derivative is given by

m ′ (t) = τ ǫ Id-e -tτǫ -1 .
This function satisfies the following lemmas :

Lemma 4.2. Let θ be an exponential random variable of mean 1. Then: The difference between jump times is denoted as

E m θ λ 0 = 1 λ 0 τ ǫ Id+λ 0 τ ǫ = A ǫ (x), (4.2
∆T c n+1 := T c n+1 -T c n , ∆T n+1 := T n+1 -T n .
The proofs of the asymptotic variance reduction make use of the following asymptotic expansions of the jump time differences of both processes : Lemma 4.4. The difference between two jump times of the control process satisfies

∆T c n+1 = θ n+1 λ 0 + ǫ θ n+1 λ 2 0 A T ǫ (X c T c n )V n + θ n+1 O(ǫ 2 ). ( 4
∆T n+1 = ∆T 0 n+1 + ǫ∆T 1 n+1 + (θ 6 n+1 + θ n+1 )O ǫ 2 + ǫnl(ǫ) , (4.4)
where

∆T 0 n+1 = θ n+1 λ 0 + b T λ 0 m ′ (∆T 0 n+1 )Z Tn = θ n+1 λ 0 + θ n+1 O(ǫ δ ) (4.5)
and, correspondingly,

∆T 1 n+1 = 1 λ 0 -b T e -∆T 0 n+1 τ -1 Z Tn b T m(∆T 0 n+1 )∇S(X Tn )V n = 1 λ 0 b T m θ n+1 λ 0 ∇S(X Tn )V n + θ n+1 O(ǫ δ ). (4.6)
It is also useful to recall that according to (2.2)-(2.13), the following hold :

|∆T n+1 | ≤ Cθ n+1 , ∆T c n+1 ≤ Cθ n+1 .
Finally, we will also need a different expansion of the jump times of the process with internal state : Lemma 4.6. When using nl 2 (ǫ) as defined by (2.11), the jump time variations of the process with internal state can be written in the following form :

∆T n+1 = θ n+1 λ 0 - b T τ ǫ λ 0 (Z Tn+1 -Z Tn ) + ǫ∆T 0 b T τ ǫ ∇S(X Tn )V n + θ 3 n+1 + θ n+1 O ǫnl 2 (ǫ) + ǫ 1+δ . (4.7)
Proof. The proof is based on Duhamel integration of the ODE (2.4) on [T n ,T n+1 ], as in the proof of Lemma 4.7 of the companion paper [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] . Following equation (4.24) in [START_REF] Rousset | Individual-based models for bacterial chemotaxis in the diffusion limit[END_REF] , we get for t ∈ [T n ,T n+1 ] :

Z t = e -(t-Tn)τǫ -1 Z Tn + ǫ Id-e -(t-Tn)τǫ -1 τ ǫ ∇S(X Tn )V n + (θ 2 n+1 + θ n+1 )O c S ǫ 2 + c F ǫ 1+δ . (4.8)
Multiplying both sides by τ ǫ , recalling that Assumption 2.4 implies τ ǫ = O(ǫ δ-1 ), and using Assumption 2.3, we get for t ∈ [T n ,T n+1 ] :

τ ǫ Z t = τ ǫ e -(t-Tn)τǫ -1 Z Tn + ǫ Id-e -(t-Tn)τǫ -1 τ ǫ 2 ∇S(X Tn )V n + (θ 2 n+1 + θ n+1 )O c S ǫ 1+δ + c F ǫ 1+γ . ( 4
Z t dt = Id-e -(Tn+1-Tn)τǫ -1 τ ǫ (Z Tn -ǫτ ǫ ∇S(X Tn )V n ) + ǫ(T n+1 -T n )τ ǫ ∇S(X Tn )V n + (θ 3 n+1 + θ 2 n+1 )O c S ǫ 2 + c F ǫ 1+δ . (4.10)
Evaluating equation (4.9) at t = T n+1 , add adding equation (4.10), we get

Tn+1 Tn Z t dt + τ ǫ Z Tn+1 =τ ǫ Z Tn + ǫ(T n+1 -T n )τ ǫ ∇S(X Tn )V n + θ 3 n+1 + θ n+1 O c S (ǫ 1+δ + ǫ 2 ) + c F (ǫ 1+δ + ǫ 1+γ
) . The estimates in Lemma 4.8 from 26 (i.e., ∆T n+1 -∆T 0 n+1 = (θ

4 n+1 + θ n+1 )O(ǫ)) yields Tn+1 Tn Z t dt + τ ǫ Z Tn+1 =τ ǫ Z Tn + ǫ∆T 0 n+1 τ ǫ ∇S(X Tn )V n + θ 3 n+1 + θ n+1 O c S (ǫ 1+δ + ǫ 2 ) + (c F + 1)ǫ 1+δ + c F ǫ 1+γ
) . Finally, plugging the last equation in the estimate (2.3) yields the result.

Analysis of variance of the coupling

The present section is devoted to the proof of the following theorem: Theorem 4.7. Assume the assumptions of Section 2.2 hold, and that kδ > 1 where k ≥ 2 is defined in (2.3). Then the difference between the process with internal state and the coupling process satisfies :

E Xt /ǫ 2 -X c t/ǫ 2 p 1/p = O(ǫ + ǫ δ + nl 2 (ǫ)), p ≥ 1, (4.11)
where nl 2 (ǫ) is defined in (2.11).

The proof of Theorem 4.7 relies on a number of steps that will be detailed by a series of lemmas. We will make use of the following random sequences, defined as the position of both processes after n jumps,

Ξ n := X Tn = X 0 + ǫ n-1 m=0 ∆T m+1 V m , n ≥ 0, (4.12)
as well as

Ξ c n := X c T c n = X 0 + ǫ n-1 m=0 ∆T c m+1 V m , n ≥ 0. (4.13)
We will also make use of the following random integers :

Definition 4.8. The random integers N ≥ 1 and N c ≥ 1 are uniquely defined by:

T N ≤ t/ǫ 2 < T N +1 , T c N c ≤ t/ǫ 2 < T c N c +1 . (4.14)
The first lemma bounds the difference at time t/ǫ 2 between the two processes Xt /ǫ 2 and X c t/ǫ 2 by expressing it in terms of differences of positions and jump times of both processes after the same random number of jumps, Ξ n and Ξ c n :

Lemma 4.9. The difference between the rescaled process with internal state X ǫ t := Xt /ǫ 2 and the rescaled coupling process X c,ǫ t := X c t/ǫ 2 satisfies:

X ǫ t -X ǫ,c t ≤ |Ξ N -Ξ c N | + ǫ|T N -T c N | + (θ N +1 + θ Nc+1 )O(ǫ) (4.15) ≤ |Ξ N c -Ξ c N c | + ǫ|T N c -T c N c | + (θ N +1 + θ Nc+1 )O(ǫ) (4.16)
Only one of the two estimates (4.15)-(4.16) is necessary in the remainder of the proof, but we detail both to highlight the symmetry.

Proof. By definition, we have

X ǫ t = Ξ N + ǫ( t/ǫ 2 -T N )V N , X c,ǫ t = Ξ c N c + ǫ( t/ǫ 2 -T c N c )V N ,
so that by Definition 4.8 of N and N c , and by realizing that t/ǫ 2 -

T c N c ≤ Cθ Nc+1 and t/ǫ 2 -T N ≤ Cθ N +1 , we get            X ǫ t -X c,ǫ t ≤ |Ξ N -Ξ c N | + |Ξ c N c -Ξ c N | (a) +(θ N +1 + θ Nc+1 )O(ǫ), X ǫ t -X c,ǫ t ≤ |Ξ N c -Ξ c N c | + |Ξ N c -Ξ N | (b) +(θ N +1 + θ Nc+1 )O(ǫ). (4.17)
To analyze (a), we consider three cases:

N = N c . Then (a) = 0. N > N c . Then (a) ǫ = N -1 n=Nc ∆T c n+1 V n ≤ T c N -T c Nc = ∆T c N c +1 + T c N -T c N c +1 ≤ T c N -T N + Cθ N c +1
, where in the last line we have used that, ∆T c N c +1 ≤ Cθ N c +1 , and, by definition,

T N ≤ t/ǫ 2 < T c N c +1 . N c > N . Then (a) ǫ = Nc-1 n=N ∆T c n+1 V n ≤ T c Nc -T c N ≤ T N +1 -T c N = ∆T N +1 + T N -T c N ≤ T N -T c N + Cθ N +1
, where in the last line we have used that by definition

T c N c ≤ t/ǫ 2 ≤ T N +1 . Variance reduced simulation of chemotaxis
This proves (4.15). By symmetry, we get:

(b) ǫ ≤ |T N c -T c N c | + C (θ N +1 + θ N c +1 ),
which yields (4.16).

In the next lemma, we estimate the supremum of the difference of the position of both processes after n jumps for n ∈ [0,n ǫ ].

Lemma 4.10. Let n ǫ ≥ 1 be a deterministic integer verifying n ǫ = O(ǫ -2 ). Then

E sup 0≤n≤nǫ |Ξ n -Ξ c n | p 1/p = O(ǫ + ǫ δ + nl(ǫ)). (4.18)
Proof. First, using the estimates of jump times (4.3)-( 4.4), we can decompose the differences of positions of both processes as follows :

Ξ n -Ξ c n = n-1 m=0 ǫ(∆T m+1 -∆T c m+1 )V m = n-1 m=0 α m+1 + d m+1 + r m+1 ,
where, by definition,

                 α m+1 := ǫ 2 Vm λ 2 0 A ǫ (Ξ m ) T -A ǫ (Ξ c m ) T V m , d m+1 := ǫ 2 Vm λ 2 0 (1 -θ m+1 )A ǫ (Ξ c m ) T V m +ǫ 2 Vm λ 2 0 λ 0 b T m θm+1 λ0 ∇S(Ξ m ) -A ǫ (Ξ m ) T V m +ǫb T m ′ (∆T 0 m+1 )Z Tm V m , r m+1 := θ 6 m+1 + θ n+1 O(ǫ 3 + ǫ 2 nl(ǫ)).
Since A ǫ is Lipschitz, we have that α m+1 ≤ Cǫ 2 |Ξ m -Ξ c m |, and we can write

|Ξ n -Ξ c n | ≤ Cǫ 2 n-1 m=0 |Ξ m -Ξ c m | + sup 0≤l≤nǫ l m=1 d m + nǫ m=1 |r m |,
for n ≤ n ǫ , which yields

|Ξ n -Ξ c n | ≤ sup 0≤l≤nǫ l m=1 d m + nǫ m=1 |r m | 1 + Cǫ 2 + ••• + (Cǫ 2 ) nǫ ,. (4.19)
Note that

1 + Cǫ 2 + ••• + (Cǫ 2 ) nǫ ≤ C.
F l = σ ((θ m ,V m-1 )|1 ≤ m ≤ l) with M 0 = 0,
so that we can apply the Burkholder-Davies-Gundy upper bound (which is a simple consequence of Doob maximal inequality here, see 1 ), for some p ≥ 1:

E sup 0≤l≤nǫ |M l | p ≤ C p n p/2-1 ǫ nǫ l=1 E(|M l -M l-1 | p ). (4.20)
In the present case, since

|Z Tn | = O(ǫ δ ): E(|d l | p ) ≤ C p ǫ p(1+δ) .
Moreover,

E nǫ m=1 |r m | p ≤ n p-1 ǫ nǫ m=1 E|r m | p ≤ C (ǫ + nl(ǫ)) p ,
and the result follows from (4.19).

In the same way, we can estimate the difference of the n ǫ -th jump time of both processes : Proof. First, using the estimates of jump times (4.3)-(4.7), we can decompose the differences of jump times of both processes as follows :

T n -T c n = n-1 m=0 ∆T m+1 -∆T c m+1 = n-1 m=0 δ m+1 + d m+1 + r m+1 ,
where, by definition,

           δ m+1 := -b T τǫ λ0 (Z Tm+1 -Z Tm ), d m+1 := ǫ∆T 0 b T τ ǫ ∇S(X Tn )V n -ǫ θn+1 λ 2 0 A T ǫ (X c T c n )V n r m+1 := θ 3 m+1 + θ n+1 ǫO(ǫ δ + nl 2 (ǫ)).
We can write for n ≤ n ǫ :

|T n -T c n | ≤ b T τ ǫ λ 0 (Z Tn -Z T0 ) + sup
E(|m m | p ) ≤ Cǫ pδ , yields E sup 0≤l≤nǫ l m=1 d m p 1/p ≤ Cǫ δ-1 .
Using:

b T τ ǫ λ 0 (Z Tn -Z T0 ) = O(ǫ 2δ-1 ) < O(ǫ δ-1 ),
and plugging in the rest term, we finally get (4.21).

We can now conclude with the proof of Theorem 4.7 :

Proof. [Proof of Theorem 4.7] We consider a given integer n ǫ ≥ 1. We decompose the difference between both processes as

X ǫ t -X c,ǫ t = X ǫ t -X c,ǫ t 1 N ≤nǫ + X ǫ t -X c,ǫ t 1 N >nǫ , (4.22)
and use (4.15) to write

X ǫ t -X c,ǫ t 1 N ≤nǫ ≤ sup 0≤n≤nǫ |Ξ n -Ξ c n | + ǫ sup 0≤n≤nǫ |T n -T c n | + (θ N +1 + θ Nc+1 )O(ǫ). ( 4
E X ǫ t -X c,ǫ t p 1 N ≤nǫ 1/p ≤ C(ǫ + ǫ δ + nl 2 (ǫ)) (4.24) + CǫE(θ p N + θ p N c ) 1/p = O(ǫ + ǫ δ + nl 2 (ǫ)), (4.25)
where in the last line we have used a technical lemma (Lemma Appendix A.1). It remains to control the probability of the event {N > n ǫ }. To this purpose, we consider t → N t ∈ N the Poisson process uniquely defined by:

Nt n=1 θ n ≤ t < Nt+1 n=1 θ n .
By definition of N and N c , there is a constant C such that N C t/ǫ 2 ≥ N and N C t/ǫ 2 ≥ N c , so that we can write:

X ǫ t -X c,ǫ t ≤ ǫC   N C t/ǫ 2 +1 n=1 θ n   1 N C t/ǫ 2 >nǫ . (4.26)
Let us choose n ǫ ≥ 4eC t/ǫ 2 , a choice that satisfies

n ǫ = O(ǫ -2 ). A second technical lemma (Lemma Appendix A.2) implies that E X ǫ t -X c,ǫ t p 1 N >nǫ ≤ ǫC p 2 -1/(Cǫ 2 ) , (4.27)
so that the contribution of the latter to the final estimate (4.11) is negligible, since it is exponentially small with respect to ǫ -2 . The proof is complete.

Discussion about the sharpness of the coupling

One may ask about the "sharpness" of the precise estimates of the coupling, given mainly by Theorem 4.7; (see also (4.18) (position shift), and (4.21) (time shift)).

Let us discuss the case where

c F = 0, δ ∈ [1/k,1/(k -1)].
In this case the dominant error term in the jump times expansion of the internal state model in Lemma 4.8 in 26 is due to the non-linearity of the turning rate: c λ ǫ kδ-1 . Since this term is due to the internal state mechanism, it can be conjectured that the difference between the density of the internal state model and gradient sensing model will be at least of this order. On the other hand, the coupling estimate in Theorem 4.7 is controlled by the same term: c λ ǫ kδ-1 , suggesting the sharpness of the latter.

Numerical illustration

In this section, we demonstrate the validity of the analysis above. For the numerical experiments, we restrict ourselves to one space dimension. For the process with internal state, we use internal dynamics that are given by the scalar cartoon model (2.5).

In all experiments, τ is chosen independently of ǫ so that δ = 1. The corresponding turning rate is given by (2.6) with λ 0 = 1. For the control process, we choose the linear turning rate (2.14), with parameters (2.22) such that the two processes have the same diffusion limit; in particular, b = 1. The physical domain is x ∈ [0,20], and we use reflecting boundary conditions, i.e. the bacterial velocity is reversed when x = 0 or x = 20. We fix a scalar bimodal chemoattractant concentration field

S(x) = α exp -β (x -ξ) 2 + exp -β (x -η) 2 , (5.1)
in which α = 5, β = 1, ξ = 7.5 and η = 12.5.

Single bacterium as a function of time. In a first experiment, we simulate a single bacterium evolving according to the fine-scale model (2.4), in which we set the parameters to ǫ = 0.2 and τ = 1. We compare this evolution to that of a bacterium that satisfies the corresponding control process (2.12). Both simulations are performed using the same random numbers starting from the initial condition X 0 = X c 0 = 8 and V 0 = V c 0 = +1; the bacterium with internal state starts from Y 0 = S(X 0 ). The time step δt = 0.1. The results are shown in figure 2. We see a very good coupling initially, which degrades over time. Note the time shift in the short time picture. In the long time picture, coupling is completely lost at some point. Variance reduced simulation of chemotaxis Expectation and variance as a function of ǫ and t. Next, we repeat the experiment using N = 10000 particles and compute the empirical mean and variance of the coupling, i.e.,

E (|Xt -X c t |) = 1 N N i X i t -X i,c t , resp., E (Xt -X c t ) 2 = 1 N N i X i t -X i,c t 2 .
As fine-scale parameters, we choose τ = 1, λ 0 = b = 1 and several values of ǫ. The chemoattractant concentration is again given as (5.1), now with α = β = 1, ξ = 7.5 and η = 12.5. The time interval for the computation is t ∈ [0,30/ǫ 2 ]. Figure 3 shows the dependence in ǫ of the coupling, by plotting the empirical mean and variance defined above as a function of ǫ for different values of the reporting time T . The results shown in figure 3 are in clear accordance with the theoretical slope predicted by the asymptotic analysis. Figure 4 shows the mean and variance of the coupling difference as a function of time. The time dependence of the variance has not been analyzed mathematically in Section 4; the specific behaviour viewed in figure 4 and is probably due to: (i) sufficiently short diffusive times; (ii) the specific double-well form of the chemoattractant potential.

Expectation and variance as a function of τ . Finally, in a last experiment, we illustrate the dependence on τ . To this end, we again simulate N = 10000 particles choosing λ 0 = 1, ǫ = 0.1 and X 0 = 7.5, for different values of τ . The results in figure 5 show that the variance quickly increases with τ until it reaches a plateau for τ > 1.

Simulation on diffusive time scales

In this section, we consider a simulation of the density of an ensemble of particles, with and without variance reduction and/or reinitialization, as described in Section 3.1. We again restrict ourselves to one space dimension, with domain x ∈ [0,20] and periodic boundary conditions. In this case, the kinetic equation corresponding to the control process reduces to the system 6. Bacterial density as a function of space at t = 50/ǫ 2 without variance reduction. Left: one realization. Right: mean over 100 realizations and 95% confidence interval. The solid line is the estimated density from a particle simulation using the process with internal state; the dashed line is estimated from a particle simulation using the control process. Both used N = 5000 particles. The dotted line is the solution of the deterministic PDE (2.20).

     ∂ t p c + + ǫ∂ x p c + = - λ c (x,+1) 2 p c + + λ c (x,-1) 2 p c - ∂ t p c --ǫ∂ x p c -= λ c (x,+1) 2 p c + - λ c (x,-1) 2 p c - . ( 6
of two PDEs, which is straightforward to simulate using finite differences.

We fix the chemoattractant concentration field as (5.1), with parameters α = 2, β = 1, ξ = 7.5 and η = 12.5. For the internal dynamics, the same model ((2.5)-(2.6), (2.14)-(2.22)) is used. The parameters are ǫ = 0.5, λ 0 = 1, τ = 1, δt = 0.1.

All simulations are performed with N = 5000 particles. The initial positions are uniformly distributed in the interval x ∈ [13,15]; the initial velocities are chosen uniformly, i.e., each particle has an equal probability of having an initial velocity of ±ǫ. The initial condition for the internal variable is chosen to be in local equilibrium, i.e., Y i 0 = S(X i 0 ). The initial positions and velocities of the control particles are chosen to be identical.

We discretize the continuum description (6.1) on a mesh with ∆x = 0.1 using a third-order upwind-biased scheme, and perform time integration using the standard fourth order Runge-Kutta method with time step δt pde = 10 -1 . The initial condition is given as

p + (x,0) = p -(x,0) = 0.25, x ∈ [13,15], 0, otherwise. (6.2)
Simulation without variance reduction. First, we simulate both stochastic processes up to time t = 50 (t = 50/ǫ 2 ) and estimate the density of each of these processes nN (x, t), resp. nc N (x, t), without variance reduction. The density is obtained via binning in a histogram, in which the grid points of the deterministic simulation are the centers of the bins. Figure 6 (left) shows the results for a single realization. We see that, given the fluctuations on the obtained density, it is impossible to conclude on differences between the two models. This observation is confirmed by computing the average density of both processes over 100 realizations. The mean densities are shown in figure 6 (right), which also reveals that the mean density of the control process is within the 95% confidence interval of the process 7. Bacterial density as a function of space at t = 50/ǫ 2 with variance reduction and reinitialization. Left: variance reduced density estimation of one realization with N = 5000 particles (solid) and deterministic solution for the control process (2.20) (dashed). Right: mean over 100 realization and 95% confidence interval (solid) and the deterministic solution for the control process (2.20) (dashed).

with internal state. Both figures also show the density that is computed using the continuum description, which coincides with the mean of the density of the control particles.

Simulation with variance reduction. Next, we compare the variance reduced estimation (3.1) with the density of the control PDE. We reinitialize the control particles after each coarse-scale step, i.e., each k steps of the particle scheme, where kδt = δt pde , (here k = 1). The results are shown in figure 7. We see that, using this reinitialization, the difference between the behaviour of the two processes is visually clear from one realization (left figure). Also, the resulting variance is such that the density of the control PDE is no longer within the 95% confidence interval of the variance reduced density estimation (right figure). We see that there is a significant difference between both models: the density corresponding to the control process is more peaked, indicating that bacteria that follow the control process are more sensitive to sudden changes in chemoattractant gradient. This difference can be interpreted from the fact that the bacteria with internal state do not adjust themselves instantaneously to their environment, but instead with a time constant τ .

Simulation with variance reduction but without reinitialization. We also compare the variance reduced estimation (3.1) with the density of the control PDE without performing any reinitialization of the control process to restore the coupling. Figure 8 (right) shows the mean estimation of the density over 200 realizations, as well as the density of the control process. As evidenced by the error bars, simulating a single realization of the process, even with variance reduction, is not able to reliably reveal this difference. This phenomenon is due to the degeneracy of the coupling on long diffusive times. This is also illustrated in figure 8 9. Bacterial density as a function of space at t = 50/ǫ 2 with variance reduction and reinitialization. Left: variance reduced density estimation of one realization with N = 5000 particles (solid) and deterministic solution for the control process (2.20) (dashed). Chemoattractant given by (5.1) with α = 5. Right: mean over 100 realization and 95% confidence interval (solid) and the deterministic solution for the control process (2.20) (dashed).

compares the variance reduced density estimation of a single realization with the density of the control PDE. Note that, as predicted by the analysis, the variance is much larger in regions where ∇S(x) is large.

Finally, we repeat the experiment with a larger chemoattractant gradient. We again choose the chemoattractant field as (5.1), with the same parameters as above, except that we now take α = 5. We also use the same discretization parameters as above. The variance reduced estimation (3.1), as well as the density of the control PDE, are shown in figure 9. We see that, although the difference between the process with internal state and the control process becomes much larger, the variance is still very well controlled by the coupling.

Conclusions and discussion

In this paper, we studied the simulation of stochastic individual-based models for chemotaxis of bacteria with internal state. We have used a coupling with a simpler, direct gradient sensing model with the same diffusion limit to obtain an "asymptotic" variance reduction. Throughout this work, we have assumed that the computation of the bacterial density using the kinetic description of the direct gradient sensing model can be performed accurately with a grid-based method. Note that in higher spatial dimensions, this may become cumbersome, and the grid-based computation should be carried out at the level of the advection-diffusion description (or any desired moment system) only, using a second level of coupling between the velocity-jump process associated with the kinetic description, and the stochastic differential equation associated with the advection-diffusion limit. This is left for future work.

The current algorithm allows to explore the differences between the fine-scale process with internal state and the simpler, coarse process using a fixed number of particles, independently of the small parameter ǫ. However, the computational cost to simulated an individual particle over diffusive time-scales becomes very large for ǫ → 0. In future work, we will therefore also study the development of truly "asymptotic preserving" schemes in the diffusion asymptotics, in the sense that the computational cost of the simulation of processes is independent of the small parameter ǫ. This will require to deal with two kinds of difficulties: (i) We will need to use an asymptotic preserving method to solve the density evolution of the control model on the grid; and (ii) We will need to extrapolate forward in time the state of the fine-scale simulation.

The first difficulty implies that, instead of solving the full kinetic equation associated with the control process, we only solve its diffusion limit. This may imply the use of a second level of variance reduction, coupling the control velocity jump process and its limiting drift-diffusion process.

The second difficulty, extrapolation in time, is related to the equation-free [START_REF] Ig Kevrekidis | Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level tasks[END_REF][START_REF] Ioannis | Equation-free multiscale computation: Algorithms and applications[END_REF] and HMM 8,9 types of methodologies; ideas of this type can be traced back to Erhenfest 10 . One approach is to use a coarse projective integration method [START_REF] Gear | coarse" integration/bifurcation analysis via microscopic simulators: micro-galerkin methods[END_REF][START_REF] Ig Kevrekidis | Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level tasks[END_REF] . One then extrapolates the bacterial density over a projective time step on diffusive time scales, after which the projected density needs to be lifted to an ensemble of individual bacteria. For such methods, besides the effect of extrapolation on the variance of the obtained results, also the effects of reconstructing the velocities and internal variables have to be systematically studied. Variance reduced simulation of chemotaxis

Step (i). The first step consists in proving the following estimate:

E θ k Nt+1 ≤ P -1 p=0 E θ k 1 θ ≥ p C t P sup p=1...P +∞ n=0 P(T n ∈ I p ) (A.2)
Indeed,

E θ k Nt+1 = +∞ n=0 E θ k n+1 1(T n ≤ t ≤ T n+1 ) ≤ +∞ n=0 P p=1 E θ k n+1 1 T n+1 -T n ≥ P -p P t|T n ∈ I p P(T n ∈ I p ).
Using the independence of θ n+1 with T n , and the assumption (A.1), it yields,

E θ k n+1 1 T n+1 -T n ≥ P -p P t T n ∈ I p ≤ E θ k 1 (θ≥ P -p CP t)
. and finally changing the sum index p to P -p + 1 yields (A.2).

Step (ii). The second step consists in the following decomposition.

.

We wish to iterate the decomposition of (a 0 ) to (a 1 ), (a 2 ), etc... For this purpose, remark that and (A.4) follows. Finally, we choose P = ⌊t + 1⌋ so that the right hand side of (A.4) does not depend on t any more. Finally we use the exponential convergence of the terms of the two series in the right hand side of (A.4) (using for instnace large deviation estimates for the second) one, to conclude that the two series indeed converge, and conclude the proof.

Lemma Appendix A.2 (Technical lemma). Let (θ n ) n≥0 be a sequence of exponential random variables of mean 1, and t → N t ∈ N the Poisson process uniquely defined by: Let n 0 ∈ N such that n 0 ≥ 4et, p 1 ≥ 1, and p 2 ≥ 1. Then there is a constant C p1,p2 independent of (n 0 ,t) such that: Step (ii). Let us denote (P t,m ,E t,m ) = (P(.|N t = m),E(.|N t = m)) the probability/expectation conditionally on the event {N t = m}. Conditionally on {N t = m}, the sequence (θ 1 ,...,θ m+1 ) has the same distribution as (tU (1) ,tU (2) -tU (1) ,...,tU (m) -tU (m-1) ,t -tU (m) + θ), where (U (1) < ••• < U (m) ) is the order statistics of m independent random variables (U 1 ,...,U m ) uniformly distributed on the interval [0,1], and θ is an independent exponential random variable of mean 1 (see Section 2.4 of [START_REF] Norris | Markov chains, volume 2 of Cambridge Series in Statistical and Probabilistic Mathematics[END_REF] for classical properties of Poisson processes). Thus denoting U (0) = 0 we can write: Step (iii). Let us compute E t,m U (n) -U (n-1)

p1p2 for n = 1,...,m + 1. First, let us recall the probability density of a couple (U (i) ,U (j) ) for 1 ≤ i < j ≤ m in an order statistics of size m:

P U (i) ∈ [u i ,u i + du i ],U (j) ∈ [u j ,u j + du j ] = m!P(U 1 < ••• < U i-1 < u i ) × P(U i ∈ [u i ,u i + du i ]) × P(u i < U i+1 < ••• < U j-1 < u j ) × P(U j ∈ [u j ,u j + du j ]) × P(u j < U j+1 < ••• < U m ),
= m! u i-1 i (i -1)! (u j -u i ) j-i-1 (j -i -1)! (1 -u j ) m-j (m -j)! du i du j . (A.7)

As a consequence,

E t,m U (n) -U (n-1) p1p2 = 0<un-1<un<1 m! u n-2 n-1 (n -2)! (u n -u n-1 ) p1p2 (1 -u n ) m-n (m -n)! du n-1 du n = m!(p 1 p 2 )! 0<un-1<un<1 u n-2 n-1 (n -2)! (u n -u n-1 ) p1p2 (p 1 p 2 )! (1 -u n ) m-n (m -n)! du n-1 du n = m!(p 1 p 2 )! (m + p 1 p 2 )! ,
where in the last line we have used (A.7) with appropriate (i,j,m) to compute the integral. In the same way, the probability density of U (i) for 1 ≤ i ≤ m in an order statistics of size m is given by:

P U (i) ∈ [u i ,u i + du i ] = m! u i-1 i (i -1)! (1 -u i ) m-i (m -i)! du i , so that E t,m U p1p2 (1)
= m!(p 1 p 2 )! Step (iv). Remarking that E m (θ p1p2 ) = (p 1 p 2 )!, and inserting the computations of Step (ii) in (A. Using the inequalities (m+1) p 2 -1 m! (m+p1p2)! ≤ 1, e

  solution Variance reduction with reinitialization Particles with internal variables Particles with gradient sensing Deterministic grid-based simulation: gradient sensing Variance reduced simulation

Fig. 1 .

 1 Fig.1. A schematic description of Algorithm 3.1. The dashed line represent the evolution of N bacteria with internal state. The dotted line represent the coupled evolution of N bacteria with gradient sensing, subject to regular reinitializations. The dashed-dotted line is computed according to a deterministic method simulating the density of the model with gradient sensing, and subject to regular reinitializations. The solid line is the variance reduced simulation of the internal state dynamics, and is computed by adding the difference between the particle computation with internal state, and the particle simulation with gradient sensing to the deterministic gradient sensing simulation. At each reinitialization step, the two simulations (deterministic and particles) of the gradient sensing dynamics are reinitialized to the values of their internal state simulation counterpart (as represented by the arrows).

Now, we can

  remark that the discrete time process (M l ) l≥0 := l m=1 d m l≥0 is a martingale with respect to the filtration

Lemma 4 . 11 .≤

 411 Let n ǫ ≥ 1 a deterministic integer. Then C p ǫ δ + nl 2 (ǫ)(4.21)

  m |. Variance reduced simulation of chemotaxis Now, we can remark that the discrete time process (M l ) l≥0 := to the filtrationF l = σ ((θ m ,V m-1 )|1 ≤ m ≤ l) with M 0 = 0.Using the Burkholder-Davies-Gundy inequality (4.20), together with:

Fig. 2 . 2 T 2 T 2 theoryFig. 3 .

 22223 Fig.2. Evolution of a single bacterium evolving according to the fine-scale process (2.4) (solid line) and the corresponding control process (2.12) (dashed) on short (left) and long (right) time-scales.

10 Fig. 4 . 2 T 2 T

 10422 Fig.4. Evolution of the empirical mean (left) and variance (right) of the difference between the fine-scale process (2.4) and the corresponding control process (2.12) as a function of t = t/ǫ 2 for different values of ǫ. The sample size is N = 10000.

Fig. 5 .

 5 Fig. 5. Empirical mean (left) and variance (right) of the difference between the fine-scale process (2.4) and the corresponding control process (2.12) as a function of τ for different values of the reporting time T . The sample size is N = 10000.

  Fig.7. Bacterial density as a function of space at t = 50/ǫ 2 with variance reduction and reinitialization. Left: variance reduced density estimation of one realization with N = 5000 particles (solid) and deterministic solution for the control process (2.20) (dashed). Right: mean over 100 realization and 95% confidence interval (solid) and the deterministic solution for the control process (2.20) (dashed).

Fig. 8 .

 8 Fig.8. Bacterial density as a function of space at t = 50/ǫ 2 with variance reduction. Left: variance reduced density estimation of one realization with N = 5000 particles (solid) and deterministic solution for the control process (2.20) (dashed). Right: mean over 200 realization and 95% confidence interval (solid) and the deterministic solution for the control process (2.20) (dashed).

  Let I = [a,b] be a finite interval of R + , T -1 = -∞. Then, +∞ n=0 P(T n ∈ I) = +∞ n=0 P(T n ∈ I,T n-1 / ∈ I) +∞ m=0 P(T n+1 ,...,T n+m ∈ I|T n ∈ I,T n-1 / ∈ I) . n-1 / ∈ I,T n ∈ I) + +∞ n=0 P(T n-1 ,T n ∈ I) = +∞ n=0 P(T n-1 / ∈ I,T n ∈ I) + +∞ n=0 P(T n ,T n+1 ∈ I)

(a m ) = +∞ n=0 PP θ 1 +P θ 1 +

 n=011 (T n ,...,T n+m ∈ I) = +∞ n=0 P(T n ∈ I)P(T n+1 ,...,T n+m ∈ I|T n ∈ I) ≤ +∞ n=0 P(T n ∈ I) P(θ 1 + ••• + θ m ≤ C(b -a)) m→+∞ -----→ 0, so that we get in the end (a 0 ) = +∞ n=0 +∞ m=0 P(T n-1 / ∈ I,T n ,...T n+m ∈ I).Factorizing P(T n ∈ I,T n-1 / ∈ I) with Bayes' formula yields (A.3).Step (iii). The estimate (A.2) yields••• + θ m ≤ t CP . (A.4)Indeed for I = [a,b], using the independence of (θ n+1 ,...,θ n+m ) from (T n ,T n-1 ):+∞ m=0 P(T n+1 ,...,T n+m ∈ I|T n ∈ I,T n-1 / ∈ I) ≤ +∞ m=0 P(θ 1 + ••• + θ m ≤ C(b -a)),and+∞ n=0 P(T n ∈ I,T n-1 / ∈ I) = P(∃n ≥ 0|T n ∈ I) ≤ 1.The two last estimates in (A.3) yields: sup p=1...P +∞ n=0 P(T n ∈ I p ) ≤ +∞ m=0 ••• + θ m ≤ t CP ,

1NN

  Nt>n0 ≤ C p1,p2 2 -n0 . Proof. November 22, 2011 12:36 WSPC/INSTRUCTION FILE variance 34 Variance reduced simulation of chemotaxisStep (i). Let us condition on the different possible values of N t > n 0 : t = m P(N t = m) t = m e -t t m m! .

5 )≤C

 5 Using the inequality (a + b) p ≤ C p a p + b b for any a,b > 0, Jensen inequality, and denoting U (m+1) = 1, there is a constant C p1,p2 such that:t -tU (m) + θ p1 + t p1 m n=1 U (n) -U (n-1) p1 p2 p1,p2 θ p2p1 + t p1p2 (m + 1)

0<u1<1 u p1p2 1 (p 1 p 2 1 =

 121 )! (1 -u 1 ) m-1 (m -1)! du m!(p 1 p 2 )! (m + p 1 p 2 )! ,and in the same way,E t,m (1 -U (m) ) p1p2 = m!(p 1 p 2 )! (m + p 1 p 2 )! .

  5)-(A.6) yields: p 2 )! + t p1p2 (m + 1) p2-1 m!(p 1 p 2 )! (m + p 1 p 2 )! e -t t m m! .

  .9) Variance reduced simulation of chemotaxis

	Moreover, integrating (4.8) on t ∈ [T n ,T n+1 ] yields Tn+1
	Tn

  -t ≤ 1 and

			+∞ m=n0+1	t m m!	≤	t n0 n 0 !	,
	it yields
	E	Nt+1 n=1	θ p1 n	p2	1 Nt>n0 ≤ C p1p2	t n0+p1p2 n 0 !	,

Rousset and Samaey

Acknowledgements

The authors thank Radek Erban, Thierry Goudon, Yannis Kevrekidis and Tony Lelièvre for interesting discussions that eventually led to this work. This work was performed during a research stay of GS at SIMPAF (INRIA -Lille). GS warmly thanks the whole SIMPAF team for its hospitality. GS is a Postdoctoral Fellow Variance reduced simulation of chemotaxis of the Research Foundation -Flanders. This work was partially supported by the Research Foundation -Flanders through Research Project G.0130.03 and by the Interuniversity Attraction Poles Programme of the Belgian Science Policy Office through grant IUAP/V/22 (GS). The scientific responsibility rests with its authors.

Appendix A. Technical lemmas, used in proof of Theorem 4.7 Lemma Appendix A.1 (Technical lemma). Let (θ n ) n≥1 a sequence of independent exponential random variables with mean 1, and (T n ) n≥1 a strictly increasing sequence of random times with T 0 = 0 such that:

Let N t a random integer such that:

Then for any integer k ≥ 0, there is a constant C k independent of t such that:

Proof. The result may seem rather intuitive, but the proof is slightly tricky. First, we consider sub-intervals of [0,t] of the form:

for some P ≥ 1 and p = 1,...,P . Variance reduced simulation of chemotaxis so that in the end, using the assumption n 0 ≥ 4et and the Stirling formula,

which yields the result.