Inferring characteristics of bacterial swimming in biofilm matrix from time-lapse confocal laser scanning microscopy

To investigate how the shape and propelling mechanism of bacteria can affect the way they navigate in a porous media such as a biofilm, we first image three bacterial swimmers –Bacillus pumilus (B. pumilus), Bacillus sphaericus (B. sphaericus), and Bacillus cereus (B. cereus) – by Transmitted Electron Microscopy (TEM) (Figure 2) to seek for potential structural and physiological differences. Important discrepancies can be observed between these Bacillus. First, they show noticeable difference in length and diameter, B. sphaericus being the longest bacteria by a factor of approximatively 1.5, and B. cereus and B. pumilus having similar size, but B. cereus showing a higher aspect ratio. Secondly, they do not have the same type of flagella: B. pumilus and B. sphaericus present several long flagella distributed over the whole surface of the membrane while B. cereus shows a unique brush-like bundle of very thin flagella, at its back tip.

Figure 2

Download asset Open asset

We then used these three species to test if these ultrastructural differences could impact their swimming behaviour in a host biofilm or in a Newtonian control fluid: could the longer body of B. sphaericus be an impediment in a crowded environment such as a biofilm or on the contrary could its larger size give it a higher strength to cross the biofilm matrix? Is the unique brush-like flagella of B. cereus an advantage or a disadvantage to swim in a Newtonian fluid or in a host biofilm?

2D+T Confocal Laser Scanning Microscopy (CLSM) of the three Bacillus swimming in a Staphylococcus aureus (S. aureus) host biofilm or in a control Newtonian buffer are acquired (see Figure 1d). Swimmers and host biofilms are imaged with different fluorescent dyes, allowing their acquisition in different color channels, and to recover in the same spatio-temporal referential the swimmer trajectories and the host biofilm density (see Materials and methods, Figure 1 and Table 1). Namely, for each species $s$ and individual swimmer $i$, we recover the initial ($T_{0,i}^s$) and final ($T_{end,i}^s$) observation times (when the swimmer goes in and out the focal plane, see Materials and methods sect. Confocal Laser Scanning Microscopy [CLSM]), and the number $T_i^s$ of time points in the trajectory. We then extract from the 2D+T images the observed position, instantaneous speed and acceleration time-series

Noting $b^s(t,x)$ the dynamic biofilm density maps obtained from the biofilm images, we also compute the local biofilm density and density gradient along trajectories

The angle ${\theta }_i^s(t)$ and the average velocity ${\overline{V}}_i^s(t)$ between two successive speed vectors are also collected (see Materials and methods sec. Post-processing of image data).

Different swimming patterns can be deciphered by qualitative observations of the trajectories $X_i^s(t)$ (Figure 3) in the biofilm and in the control Newtonian buffer, and run-and-tumble swimming patterns are quantified with ${\theta }_i^s(t)$ and ${\overline{V}}_i^s(t)$ (Figure 4). B. sphaericus has a similar run-and-reverse behaviour in the biofilm and the control buffer with trajectories divided between back and forth paths around the starting point and long runs, the biofilm strongly impairing its speed and increasing the number of reverse events. By contrast, B. pumilus clearly switches its swimming behaviour in the biofilm, from quasi-straight runs in the Newtonian buffer to a pronounced run-and-reverse behaviour in the biofilm with decreased speeds and chaotic trajectories. On the contrary, B. cereus swimmers manage to conserve comparable trajectories and distributions of swimming speed and direction in the biofilm compared to control. Interestingly, the number of reverse events is even reduced in the host biofilm for B. cereus.

Figure 3

Download asset Open asset

Figure 4

Download asset Open asset

For further quantitative analysis, trajectory descriptors are defined. We first investigate the distribution of the population-wide average acceleration and velocity norms ${\displaystyle \frac{1}{T_i^s-2}{\sum }_t\Vert A_i^s(t)\Vert }$ and ${\displaystyle \frac{1}{T_i^s-1}{\sum }_t\Vert V_i^s(t)\Vert }$, where $\parallel \cdot \parallel$ denotes the Euclidian norm. We also quantify the swimming kinematics by computing the travelled distance $dis{t}_i^s$ along the path and the total displacement $dis{p}_i^s$, that is the distance between the initial and final trajectory points, with

We finally compute the total biofilm area visited by a swimmer along its path (see Figure 1b). The same descriptors are computed in the control Newtonian buffer.

The three species present contrasted distributions for these descriptors (Figure 5). B. sphaericus has the smallest mean ($||A||=0.58$ and $||V||=0.70$) and median ($\parallel A\parallel =0.50$ and $\parallel V\parallel =0.53$) values of acceleration and speed, while B. pumilus has the widest distributions (difference between 95% and 5% centiles of 2.76 for $\parallel A\parallel$ and 2.45 for $\parallel V\parallel$ compared to 1.00, 1.51 and 1.90, 1.49 for B. sphaericus and B. cereus respectively). B. cereus for its part shows the highest accelerations, indicating larger changes in swimming velocities, but median and mean speeds comparable to B. pumilus (Figure 5, $\parallel A\parallel$ and $\parallel V\parallel$ panels). We also note that B. sphaericus and to a lower extent B. pumilus trajectories have a significant amount of null or small average speeds, while B. cereus trajectories have practically no zero velocity, consistently with the qualitative analysis (Figure 5, $\parallel V\parallel$ panels). Small velocities episodes of B. sphaericus and B. pumilus could occur during their back-and-forth trajectories, which produce small displacements and pull the displacement distribution towards lower values than B. cereus (Figure 5, Disp panel). B. pumilus displacement is intermediary. Conversely, back-and-forth trajectories can produce large swimming distances for B. sphaericus and B. pumilus (mean adimensioned value of 32.2 and 43.2 respectively) so that B. sphaericus has a distance distribution comparable to B. cereus (mean adimensioned value of 29.6, Figure 5, Dist panel), but lower than B. pumilus. Observing conjointly displacement and distance (Figure 5, lower-right panel) provides consistent insights: B. sphaericus shows a large variability of small displacement trajectories, from small to large distances, while B. cereus trajectory displacement seems to vary almost linearly with the distance at least for the points inside the isoline 50%. B. pumilus has again an intermediary distribution, with a large range of displacement-distance couples. The distributions of visited areas of B. pumilus and B. cereus are almost identical, and higher than B. sphaericus one. Compared to the control buffer, all descriptors are reduced in the biofilm. Consistently with previous observations, the displacement ($disp$) is strongly reduced for B. pumilus, and less impacted for B. sphaericus and B. cereus. These observations must be related to the behavioural switch for B. pumilus and to the identical swimming patterns for the two other Bacilii in the biofilm compared to the control fluid.

Figure 5

Download asset Open asset

All together, this data depict (1) a long-range species, B. cereus, which moves efficiently in the biofilm during long, relatively straight, rapid runs, almost identically as in a Newtonian fluid (2) a short-range species, B. sphaericus, that moves mainly locally in small areas in the biofilm and in the control buffer with lower accelerations and speeds except few exceptions (only 6% of its trajectories induced a displacement higher than $10\mu m$ compared to 28% for B. cereus and 26% for B. pumilus) and (3) a medium-range species, B. pumilus, with a large diversity of rapid trajectories, from small to large displacement, and a behavioural change from straight runs in a Newtonian fluid to frequent run-and-reverse events in the biofilm. These kinematics discrepancies for B. pumilus and B. cereus allow them however to cover identical visited areas.

Though, these global descriptors do not inform about potential adaptations of the swimmers to the biofilm matrix. We first check if swimmer velocities are directly linked to the local biofilm density, and if the swimmers adapt their trajectory according to density gradients by plotting the points $(\parallel \nabla b(t,X_i^s(t))\parallel ,\parallel A_i^s(t)\parallel )$ and $(b(t,X_i^s(t)),\parallel V_i^s(t)\parallel )$ (Figure 5, lower panel). Clear differences between the three species can be deciphered. First, the three Bacillus do not have the same distribution of visited biofilm density and gradient. B. pumilus swimmers visit denser biofilm with higher variations than the other species while B. sphaericus and B. cereus stay in less dense and smoother areas, the quantile 0.5 of these species being circumscribed in low gradient and low density values. Next, B. cereus has a wider distribution of accelerations, specially for small-density gradients, compared to B. pumilus and B. sphaericus. This could indicate that when the biofilm is smooth, B. cereus samples its acceleration in a large distribution of possible values. Finally, we observe that the speed distribution rapidly drops for increasing biofilm densities for B. sphaericus and B. cereus, while the decrease is much smoother for B. pumilus. These observations provide additional insights in the species swimming characteristics: B. pumilus swimmers seem to be less inconvenienced by the host biofilm density than the other species, while B. cereus and B. sphaericus bacteria appear to be particularly impacted by higher densities and to favor low densities where it can efficiently move. Though, B. sphaericus has lower motile capabilities than B. cereus when the biofilm is not dense.

This descriptive analysis does not allow to clearly identify potential mechanisms by which the swimmers adapt their swim to the biofilm structure or to simulate new species-dependant trajectories. We then build a swimming model based on a Langevin-like equation on the acceleration that involves several swimming behaviours modelling the swimmer adaptation to the biofilm. Furthermore, after inference, new synthetic data can be produced by predicting swimmer random walks sharing characteristics comparable to the original data.

We consider bacterial swimmers as Lagrangian particles and we model the different forces involved in the update of their velocity $\mathbf{v}$. We assume that the swimmer motion can be modelled by a stochastic process with a deterministic drift (Figure 1c):

where the right hand side is composed of two deterministic terms in addition to a gaussian noise, each weighted by the parameters $\gamma$, $\beta$ and $ϵ$.

The first term implements the biological observation (Figure 5, lower central panel) that the bacterial swimmers adapt their velocity to the biofilm density. This term can be interpreted as a speed selection term that pulls the instantaneous speed of the swimmer towards a prescribed target velocity $\alpha (b)$ that depends on the host biofilm density $b$. The weight $\gamma$ can be interpreted as a penalization coefficient. In such a formalism, the difference between the swimmer and the prescribed speed is divided by a relaxation time $\tau$ to be homogeneous to an acceleration. Hence, $\gamma$ is proportionally inverse to $\tau$, $\gamma \sim \frac{1}{\tau }$. As a first-order approximation of the speed drop observed in Figure 5 for increasing $b$, the target speed $\alpha (b)$ is modeled as a linear variation between v₀ and v₁, where v₀ is the swimmer characteristic speed in the lowest density regions, where $b=0$, and v₁ in the highest density zones where $b=1$:

The second term updates the velocity direction according to the local gradient of the biofilm density $\nabla b$. The sign of $\beta$ indicates if the swimmer is inclined to go up (negative $\beta$) or down (positive $\beta$) the host biofilm gradient, while the weight magnitude indicate the influence of this mechanism in the swimmer kinematics. We note that this term does not depend on the gradient magnitude but only on the gradient direction: this reflects the implicit assumption that the bacteria are able to sense density variations to find favorable directions, but that the biological sensors are not sensitive enough to evaluate the variation magnitudes.

The third term is a stochastic two-dimensional diffusive process that models the dispersion around the deterministic drift modelled by the two first terms. We define

The term η can also be interpreted as a model of the modelling errors, tuned by the term $ϵ$. Equation 1 is supplemented by an initial condition by swimmer. For vanishing $\parallel v\parallel$ or $\parallel \nabla b\parallel$ leading to an indetermination, the corresponding term in the equation is turned off.

Equation 1 links the observed biofilm density and the swimmer trajectories trough mechanistic swimming behaviours. The model fitting can be seen as an ANOVA-like integrative statistical analysis of the image data. It decomposes the observed acceleration variance between mechanistic processes describing different swimming traits in order to decipher their respective influence on the swimmer trajectories while integrating heterogeneous data (density maps $b$ and trajectories kinematics).

We can define characteristic speed and acceleration $V^{*}$ and $A^{*}$ in order to set a dimensionless version of Equation 1

where ${\gamma }^{\prime }=\frac{\gamma V^{*}}{A^{*}}$, $v_0^{\prime }=\frac{v_0}{V^{*}}$, $v_1^{\prime }=\frac{v_1}{V^{*}}$, ${\beta }^{\prime }=\frac{\beta }{A^{*}}$, ${\displaystyle {\eta }^{\mathrm{\prime }}\sim \mathcal{N}(0,{ϵ}^{\mathrm{\prime }})}$ and ${ϵ}^{\prime }=\frac{ϵ}{A^{*2}}$.

This dimensionless version will strongly improve the inference process and will allow an analysis of the relative contribution of the different terms in the kinematics. An extended numerical exploration of this model is performed in Appendix 2 Sec. Numerical exploration on mock biofilm images to illustrate the impact of the different parameters on the trajectories, showing in particular the interplay between $\gamma$ and $ϵ$: counter-intuitively, straight lines are induced when the stochastic part $ϵ$ is high compared to the speed selection parameter $\gamma$ (see also Appendix 2).

For each bacterial swimmer population, we now seek to infer with a Bayesian method population-wide model parameters governing the swimming model of a given species from microscope observations.

Inference model setting

Equation (2) is re-written as a state equation on the acceleration for the bacterial strain $s$ and the swimmer $i$

(3) $\begin{array}{cc}\hfill A_i^s(t)& ={\gamma }^s(v_0^s+b(t,X_i^s(t))(v_1^s-v_0^s)-\parallel V_i^s(t)\parallel )\frac{V_i^s(t)}{\parallel V_i^s(t)\parallel }+{\beta }^s\frac{\nabla b(t,X_i^s(t))}{\parallel \nabla b(t,X_i^s(t))\parallel }+{\eta }^s\hfill \end{array}$

(4) $\begin{array}{cc}& :=f_A({\theta }^s,b(t,X_i^s(t)),V_i^s(t),X_i^s(t))+{\eta }^s\hfill \end{array}$

where

${\theta }^s:=({\gamma }^s,v_0^s,v_1^s,{\beta }^s)$

are species-dependant equation parameters. The function $f_A$ can be seen as the deterministic drift of the random walk, gathering all the mechanisms included in the model. The inter-individual variability of the swimmers of a same species comes from the swimmer-dependent initial condition, the resulting biofilm matrix they encounter during their run, and the stochastic term.

Inferring the parameters ${\theta }^s$ can then be stated in a Bayesian framework as solving the non linear regression problem

(5) $\begin{array}{r}{A}_i^s(t)\sim \mathcal{N}\left(f_A\left({\theta }^s|b(t,X_i^s(t)),V_i^s(t),X_i^s(t)\right),{ϵ}^s\right)\end{array}$

from the data $b(t,X)$, $X_i^s(t)$, $V_i^s(t)$ and $A_i^s(t)$, with truncated normal prior distributions

(6) $\begin{array}{r}{\theta }^s\sim \mathcal{N}(0,1),{ϵ}^s\sim \mathcal{N}(0,1),\end{array}$

and additional constrains on the parameters

${\gamma }^s\ge 0,v_0^s\ge 0,v_1^s\ge 0,{ϵ}^s\ge 0.$

We note that Equation (5) can be seen as a likelihood equation of the parameter ${\theta }^s$ knowing $A_i^s(t),b(t),V_i^s(t)$ and $X_i^s(t)$. The parameter ${ϵ}^s$ can now be seen as a corrector of both modelling errors in the deterministic drift and observation errors between the observed and the true instantaneous acceleration. Alternative settings where these uncertainties sources are separated and a true state for position and acceleration is inferred can be defined (see Annex Various inference models). The inference problem is implemented in the Bayesian HMC solver Stan (Stan Development Team, 2018) using its python interface pystan (Riddell et al., 2021). Inference accuracy is thoroughly assessed on synthetic data (see Appendix 1 Assessment of the inference with synthetic data and Figure 6).

Figure 6

Download asset Open asset

Analysis of the confocal microscopy dataset

We now solve the inference problem (5)-(6) on the confocal microscopy dataset to identify population-wide swimming model parameters in order to decompose the swimmer kinematics in three mechanisms: biofilm-related speed selection, density-induced direction changes and random walk. The inference process is assessed by comparing the descriptors obtained on trajectories predicted by the fitted model (Figure 7a) with descriptors of real trajectories (Figure 5). The mean values of acceleration and speeds are accurately predicted for the three species ( Figure 7a panels $\parallel A\parallel$ and $\parallel V\parallel$, dashed lines). Relative positions of distance, displacement and visited area mean values are also correctly simulated (Figure 5 and Figure 7a, upper panel). B. sphaericus presents the lowest predicted accelerations and speeds while B. pumilus has the widest speed and acceleration distributions and B. cereus shows the highest accelerations, consistently with the data. The visited area and the distances are slightly over estimated, but the relative position and the shape of the distributions are conserved. The amount of null velocities for B. sphaericus is under estimated by the fitted model and not rendered for B. pumilus. The distance distributions of the three species are accurately predicted by the fitted model. When displaying conjointly the distance and the displacement (Figure 7a, right lower panel), the distribution of B. sphaericus is correctly predicted by the simulations, but B. cereus and B. pumilus displacements are underestimated. Some qualitative features can be recovered, such as the higher distribution of distance-distribution couples for B. cereus or higher displacement for B. cereus compared to B. sphaericus.

Figure 7

Download asset Open asset

Descriptors of swimming adaptations to the host biofilm are also correctly preserved for the main part (Figure 5 and Figure 7 a, lower panel). B. pumilus is the species that crosses the highest biofilm densities in the fitted model simulations, showing the highest speeds in this crowded areas, and that visits the most frequently areas with high density gradients, consistently with the data. As in the confocal images, the simulated B. sphaericus and B. cereus favor smoother zones of the biofilm with lower biofilm densities. The B. cereus fitted model correctly renders the highest acceleration variance observed in the data for low biofilm gradients, while B. sphaericus speed and acceleration variance is the lowest for all ranges of biofilm densities and gradients, both in the data and in the fitted model predictions. The drop of speeds and accelerations for increasing biofilm densities and gradients is well predicted for B. pumilus, but is smoother in the simulation compared to the data for B. sphaericus and B. cereus. In particular, the sharp drop of speeds for $b\simeq 0.25$ observed in the data for B. cereus and B. sphaericus is underestimated by the fitted model.

All together, the model reproduces very accurately the mean values of acceleration, speed and visited area, renders relative positions and the main characteristics of distributions for distance, displacement and interactions with the host biofilm matrix, but produces less variable outputs than observed in the data, meaning that the model is less accurate in the distribution tails. The main features of the swimmer adaptation to the underlying biofilm are however correctly predicted by the model.

To further inform the fitted model accuracy, the coefficient of determination $R_{det}^2$ of the deterministic components $f_A({\theta }^s,b(t),V_i^s,X_i^s(t))$ of Equation 4 is computed (Table 2), in order to quantify the goodness of fit of the friction and gradient terms of (Equation 2) that represent interactions with the biofilm. These results highlight that B. cereus bacteria do present an important stochastic part in the accelerations, while the B. pumilus species is the best represented by our deterministic modelling.

Table 2

data	N	$A_{ref}$	$V_{ref}$	$\sigma (A)$	$R_{det}^2[\%]$	${ϵ}^2$
B.pumilus	33,916	81.08	7.89	0.87	58.80	0.36
B.sphaericus	20,152	44.93	4.74	0.58	48.50	0.30
B.cereus	23,160	108.92	7.03	0.63	32.72	0.42

The three species present very different inferred parameter values (Figure 7 b and Table 3), showing that the model inference captures contrasted swimming characteristics of these Bacillus. Due to the mechanistic terms introduced in Equation 1, these differences can be interpreted in term of speed and direction adaptations to the host biofilm. First, B. pumilus shows the highest v₀ value, and the highest amplitude between v₀ and v₁, inducing a higher ability for B. pumilus to swim fast in low density biofilm zones and strong deceleration in crowded area. In comparison, B. sphaericus presents the smallest amplitude between v₀ and v₁ showing a poor adaptation to biofilm density. B. cereus has the highest $\gamma$ value, showing a reduced relaxation time toward the density dependant speed: in other words, B. cereus is able to adapt its swimming speed more rapidly than the other species when the biofilm density varies. B. cereus swimmers are also better able to change their swimming direction in function of the biofilm variations they encounter along their way, their $\beta$ distribution being markedly higher than the other species which have very low $\beta$. Finally, the stochastic parameter $ϵ$ is also contrasted, from a low distribution for B. sphaericus to high values for B. cereus. All together, the inference complete the observations made in Figure 5: B. pumilus poorly adapts its swimming direction to the host biofilm (low $\beta$) but has a wide range of possible speeds when the biofilm density varies (high v₀, low v₁), that it can reach quite rapidly (intermediary $\gamma$) with intermediary stochastic correction ($ϵ$). In contrast, B. cereus reaches lower speed values (intermediary v₀, low v₁) but is more agile to adapt its swimming to its environment by changing rapidly its speed when the biofilm density is more favorable (highest $\gamma$) and adapting its swimming direction to biofilm variations, with higher stochastic variability (large $ϵ$). Finally, B. sphaericus is the less flexible of the three bacteria: less fast (smallest difference between v₀ and v₁), they are also less responsive to biofilm variations (small $\gamma$ and $\beta$) with low random perturbations (small $ϵ$).

Table 3

species	param	mean	std	confidence interval [2.5%–97.5%]	neff	Rhat
B. pumilus	γ	0.77	3.95×10–3	[0.77–0.77]	4,507	1
	v0	0.14	8.67×10–3	[0.12–0.16]	3,879	1
	v1	1.69×10–3	1.69×10–3	[5.18×10–5−6.26×10–3]	4,821	1
	β	9.84×10–3	5.07×10–3	[1.45×10–5−2.07×10–2]	5,223	1
	ε	0.62	2.48×10–3	[0.61–0.62]	5,307	1
B. sphaericus	γ	0.61	4.53×10–3	[0.60–0.62]	4,965	1
	v0	2.75×10–4	2.75×10–4	[4.91×10–6−1.01×10–3]	4,019	1
	v1	4.84×10–3	4.77×10–3	[9.39×10–5−1.45×10–2]	5,001	1
	β	4.25×10–3	3.33×10–3	[−2.18×10–3−1.15×10–2]	4,668	1
	ε	0.32	1.55×10–3	[0.31–0.32]	5,943	1
B. cereus	γ	0.83	1.11×10–2	[0.80–0.86]	2,700	1
	v0	6.44×10–2	1.07×10–2	[3.22×10–2−9.66×10–2]	2,510	1
	v1	6.65×10–3	6.33×10–3	[1.50×10–4−2.15×10–2]	4,061	1
	β	2.78×10–2	9.04×10–3	[1.39×10–2−5.56×10–2]	4,230	1
	ε	0.9	4.17×10–3	[0.89–0.92]	4,852	1

Finally, after inference, the impact of each term in the overall acceleration data can be quantified and analyzed by displaying its relative contribution in a ternary plot (Appendix 2—figure 6). This relative contribution can be measured thanks to the swimming model which integrates these different mechanisms in the same inference problem. The direction selection is the least influential mechanism for the three species, with a slightly higher impact for B. cereus (50% and 95% isolines slightly shifted towards $A(\nabla b)$ in Appendix 2—figure 6a). When zooming in, the three Bacillus show differences in the balance between speed selection and the random term (Appendix 2—figure 6b): while B. pumilus is slightly more influenced by the friction term than by stochasticity, these mechanisms are perfectly balanced in B. sphaericus accelerations, while B. cereus is more influenced by the random term.

When analyzing microbial swimming trajectories, two general strategies can be found in the literature. The first one aims at designing statistical tests quantifying similarities with or deviations from typical motion of interest such as diffusion (Patteson et al., 2015). Another strategy consists in providing a generative model of the data, analyzing it (Chepizhko and Peruani, 2013; Chepizhko et al., 2013) and comparing model outputs with real data (Koorehdavoudi et al., 2017; Jabbarzadeh et al., 2014), possibly after inference. The model that is studied in this paper belong to the second category: the model includes deterministic mechanisms describing interactions with the host biofilm, together with a random correction counterbalancing the modelling errors. The parameter inference allows to interpret the data variance relatively to speed or direction adaptations to the host biofilm versus residual effects gathered in the stochastic term. This method is comparable to ANOVA-like multivariate analysis: the parametric phenomelogical mappings between explicative co-variables and a swimming behaviour (for example the function defining speed selection from biofilm density) are gathered in the same inference problem, enabling to decompose acceleration variability between the different swimming behaviours. This integrative method allows for multi-data integration and co-analysis. Furthermore, the fitted model allows to simulate typical swimming trajectories of a given species.

In this study, we do not aim to recover ’true’ swimmer trajectories (a.e. the blue trajectory in Appendix 2—figure 4), that is identifying through smoothing techniques an approximation of the specific realization of the stochastic modeling and observation errors that lead to a given ’observed’ trajectory. Rather, the goal is to identify common characteristics shared by a population of trajectories by inferring the ‘population-wide’ parameters (the parameters $\alpha$, $\beta$, v₀, v₁, $\gamma$, and $ϵ$) that best explain the whole set of observed accelerations in a same population of swimmers. For this reason, we did not introduce swimmer-specific terms nor individual noise: they would have increased the model accuracy, but to the price of a blurrier characterization of the species specificities.

This choice determined our inference framework. Despite several alternative options for recovering hidden states, in particular SSM (space state models) which are common in spatial ecology (Auger‐Méthé et al., 2021), the Bayesian method we opted for is a simpler non-linear regression problem that proved to be sufficient to recover macroscopic features of swimmer trajectories and species stratification. We discuss in Appendix 3 Various inference models the different options that were tested and present in Materials and methods Sec. Inference the method for noise model selection. Among other interesting features, the Bayesian method provides confidence intervals on the final parameter estimation, and on the resulting trajectories as in Figure 6a.

The deterministic terms of the model explain only half of the variance (Table 2). A major part of the underlying mechanisms is not correctly described by our model which is a common feature since it is a phenomenological model which only considers interactions with the underlying biofilm at a macroscopic level, without taking into account nanoscale physical mechanisms. A more detailed description of the underlying physics could have been designed as in Martinez et al., 2014, but it would have made more complex the analysis of the interactions between the host biofilm and the swimmer trajectories and the extraction of species-specific patterns. However, we note that our model correctly renders observations made through macroscopic trajectory descriptors, even though the inference process has not been made based on these observables. Furthermore, several repetitions of the same models with different samples of the stochastic terms give very similar values for the trajectory descriptors (see Appendix 2—figure 5 and section Influence of inference and stochastic terms on the trajectory descriptors), showing that these descriptors are robust to stochastic perturbations. Hence, the model (2) can be used to produce synthetic data sharing the same global characteristics than the original ones specifically taking into accounts interactions between the swimmers and the host biofilm. Furthermore, these predictions also reproduce the species stratification observed in the original data using the global descriptors.

The direction selection term of the equation driven by $\beta$ has little impact in the swimmer model fitted on real data. However, the parameter $\beta$ can have a sensible impact on the kinematics as shown in the sensitivity analysis, and on the trajectories in mock biofilms (Appendix 2—figure 1). This could indicate that direction selection based on biofilm gradients is marginally effective in real-life swimming trajectories in a biofilm matrix. On the contrary, the speed selection term is more effective for the three Bacillus, showing that these micro-swimmer are able to adapt their swimming velocity to the biofilm density faced during their run. This term acts as an inertial term which enhances the stochastic term to provide direction and velocity changes.

The model has been used to decipher different adaptation strategies to the host biofilm of the three species during their swim. It confirms that B. sphaericus are the less motile bacteria in the biofilm, with reduced speeds and adaptation capabilities as indicated by the smallest model parameter values and a stereotypic run-and-reverse behaviour inside or outside the biofilm. B. pumilus on the contrary drastically changes its swimming behaviour in the biofilm compared to the Newtonian control buffer, which is reflected in the model by a high amplitude between v₀ and v₁ and a high $\gamma$ that indicates a rapid adaptation for varying biofilm densities. B. cereus shows the highest adaptation ability to the biofilm matrix, with the highest $\gamma$ and $\beta$ reflecting biofilm-induced speed and direction changes. Furthermore, the high stochastic effects (highest ∈) higher than the speed selection term tuned by $\gamma$ (see Appendix 2—figure 6) allows this swimmer to conserve straight runs in the biofilm (see Appendix 2 Sec. Friction and random term in Langevin equations.) in the same way than in the control Newtonian fluid.

This characterization methodology could be used to drive species selection for improved biofilm control. Furthermore, the model can be used to predict new trajectories and the resulting biofilm vascularization, in a similar framework as in Houry et al., 2012. Coupled with a model of biocide diffusion, these simulations could be used to test numerically the efficiency of mono- or multi-species swimmer pre-treatment to improve the removal of the host biofilm.

Characterization of flagellated bacteria motility in polymeric solutions is a very active research area (Martinez et al., 2014; Patteson et al., 2015; Zöttl and Yeomans, 2019; Qu and Breuer, 2020; Qu et al., 2018). Speed and direction variations have been measured for various polymeric fluids with different visco-elastic properties. For the model bacteria E. coli in polymeric solutions, enhanced viscosity decreases tumbling while increased elasticity speeds up the swimmers (Patteson et al., 2015; Zöttl and Yeomans, 2019). In our experiments on the contrary, we observed decreased speeds and strong enhancement of reverse events for the flagellated B. sphaericus and B. pumilus in the biofilm compared to the Newtonian control buffer. However, the experimental set-up shows strong differences: the complex rheology of S. aureus biofilms may strongly differ from polymeric fluids even if under certain condition they can be considered as visco-elastic fluids (Gloag et al., 2020), impacting differently the swimmer behaviours. Furthermore, the physiology of the motor cell in the Gram-positive Bacillus differs from the one of the Gram-negative E. coli (Terahara et al., 2020; Szurmant and Ordal, 2004; Subramanian and Kearns, 2019). Finally, the particular brush-like flagella bundle of B. cereus may allow this species to conserve the same swimming in Newtonian and crowded environments, by adapting its swimming speed to the local density and otherwise randomly selecting swimming directions across the host biofilm. To generalize this approach to other contexts, this study should be reproduced for other swimmers and other host biofilms, together with polymeric fluids and porous media, including biochemical interactions.

Infiltration of S. aureus biofilms by bacilli swimmers were prepared in 96-well microplates. Submerged biofilms were grown on the surface of polystyrene 96-well microtiter plates with a μ clear base (Greiner Bio-one, France) enabling high-resolution fluorescence imaging (Bridier et al., 2010). 200 μL of an overnight S. aureus RN4220 pALC2084 expressing GFP (Malone et al., 2009) cultured in TSB (adjusted to an OD 600 $\mathrm{nm}$ of 0.02) were added in each well. The microtiter plate was then incubated at 30°C for 60 $\min$ to allow the bacteria to adhere to the bottom of the wells. Wells were then rinsed with TSB to eliminate non-adherent bacteria and refilled with 200 μL of sterile TSB prior incubation at 30 $\mathrm{celsius}$ for 24 h. In parallel, B. sphaericus 9 A12, B. pumilus 3 F3 and B. cereus 10B3 were cultivated overnight planktonically in TSB at 30 °C. Overnight cultures were diluted 10 times and labelled in red with 5 μM of SYTO 61 (Molecular probes, France). After 5 min of contact, 50 μL of labelled fluorescent swimmers suspension were added immediately on the top of the S. aureus biofilm. All microscopic observations were collected within the following 30 min to avoid interference of the dyes with bacterial motility. Three replicates were conducted. The same protocol has been repeated without the host biofilm (control experiments): the swimmers are added to the buffer only which is a Newtonian fluid.

The 96 well microtiter plate containing 24 hr S. aureus biofilm and recently added bacilli swimmers were mounted on the motorized stage of a Leica SP8 AOBS inverter confocal laser scanning microscope (CLSM, LEICA Microsystems, Germany) at the MIMA2 platform (https://www6.jouy.inra.fr/mima2_eng/). Temperature was maintained at 30 $\mathrm{celsius}$ during all experiments. 2D+T acquisitions were performed with the following parameters: images of 147.62 × 147.62 $\mu \mathrm{m}$ were acquired at 8000 $\mathrm{Hz}$ using a 63×/1.2 N.A. To detect GFP, an argon laser at 488 $\mathrm{nm}$ set at 10% of the maximal intensity was used, and the emitted fluorescence was collected in the range 495–550 $\mathrm{nm}$ using hybrid detectors (HyD LEICA Microsystems, Germany). To detect the red fluorescence of SYTO61, a 633 $\mathrm{nm}$ helium-neon laser set at 25% and 2% of the maximal intensity was used, and fluorescence was collected in the range 650–750 $\mathrm{nm}$ using hybrid detectors. Images were collected during 30 s (see Table 1 for sampling period).

Bacterial swimmers navigate within a three-dimensional biofilm matrix and confocal microscope refreshment time is not small enough to allow 3D+T images. To limit 3D trajectories, a focal plane near the well edge has been selected, where the well wall physically constrains the swimmer trajectories in one direction, which select longer trajectories in the 2D plane that can be tracked in time. Therefore, experimental data are composed of two-dimensional trajectories captured between the swimmer arrival and departure times in the focal plane, and the associated 2D+T biofilm density images that change over time due to swimmer action.

To check that the host biofilm structure is identical near the well’s edge compared to other 2D slices, we took 4 replicates of S. aureus biofilms that were imaged in 3D using a stack of 6 horizontal images, starting from $z=0$ near the well’s edge, to $z=6\mathrm{\Delta }z$, at the interface between the biofilm and the bulk solution. To study the between and within biofilm density variability in the horizontal images, we subsampled them with a regular Cartesian 4 × 4 grid, resulting in a 4 × 6 x(4 × 4)=384 2D images database supplemented by metadata (stack, $z$ and $x-y$ coordinate of the subsample), before computing a clustered pairwise correlation similarity matrix and a permanova.

Materials were directly adsorbed onto a carbon film membrane on a 300-mesh copper grid, stained with 1% uranyl acetate, dissolved in distilled water, and dried at room temperature. Grids were examined with Hitachi HT7700 electron microscope operated at 80 $\mathrm{kV}$ (Elexience – France), and images were acquired with a charge-coupled device camera (AMT).

See Figure 1 for a sketch of the datastream from microscope raw images to model inputs and Appendix 1—figure 1 for data visualization at each step of the post-processing pipeline.

Swimmer tracking has been applied on the red channel of the raw temporal stacks with IMARIS software (Oxford Instruments) using the tracking function after automated spots detection to get position time-series for each swimmer. Time-series with less than 8 time steps were filtered out.

Then, swimmer speed and acceleration time-series were computed from their position by finite-difference approximations and trajectory descriptors were extracted. The RGB green channel corresponding to the biofilm density temporal images were converted into grayscale and rescalled between 0 and 1 (linear scalling).

Trajectory descriptors are defined as follows. The mean acceleration and speed values, distance and displacement are computed with ${\parallel A\parallel }_i^s=\frac{1}{T_i^s-2}{\sum }_t\parallel A_i^s(t)\parallel$, ${\parallel V\parallel }_i^s=\frac{1}{T_i^s-1}{\sum }_t\parallel V_i^s(t)\parallel$, $dis{t}_i^s=\mathrm{\Delta }t{\sum }_{T_{0,i}^s}^{T_{end,i}^s-\mathrm{\Delta }t}\parallel V_i^s(t)\parallel$ and $dis{p}_i^s=\parallel X(T_{end,i}^s)-X(T_{0,i}^s)\parallel$. To compute the visited area, each trajectory piece was subsampled by computing $X_i^s(t_k)=\frac{k}{n_s}{X}_i^s(t)+(1-\frac{k}{n_s})X_i^s(t+\mathrm{\Delta }t)$ for $k=0,n_s$, with $n_s=10$ and the pixels included in the ball $B(X_i^s(t_k),r)$ with radius $r=2$ were labeled. The total area of the labelled pixels is defined as the visited area of the swimmer $i$ of species $s$.

To assess run-and-tumble behaviour, the angle ${\theta }_i^s(t)$ and the mean velocity ${\overline{V}}_i^s(t)$ between two consecutive speed vectors are defined with ${\displaystyle {\theta }_i^s(t)}$ = ${\displaystyle \arccos ((V_i^s(t)\cdot V_i^s(t-\mathrm{\Delta }t))/(\Vert V_i^s(t)\Vert \Vert V_i^s(t-\mathrm{\Delta }t)\Vert ))}$ and ${\displaystyle {\overline{V}}_i^s(t)=(\Vert V_i^s(t)\Vert +\Vert V_i^s(t-\mathrm{\Delta }t)\Vert )/2}$, for ${\displaystyle t\in (T_{0,i}^s+\mathrm{\Delta }t,T_{end,i}^s)}$.

Post-processed data are available at https://forgemia.inra.fr/bioswimmers/swim-infer/SwimmerData.

Time integration of equations (2) has been solved with an explicit Euler scheme regarding positions ${\mathbf{x}}_{i,t}^s$ and velocities ${\mathbf{v}}_{i,t}^s$ of the swimmer $i$ of species $s$ at time $t$:

where ${\mathbf{dv}}_{i,t}^s$ is given by Equation 2, and depends on ${\theta }^s$, $V_{i,t}^s$, $x_{i,t}^s$, $b(t,x_{i,t}^s)$ and $\nabla b(t,x_{i,t}^s)$. In practice, the biofilm density and gradient maps $b$ and ${\displaystyle \mathrm{\nabla }b}$ are discretized with a Cartesian grid corresponding to the image pixels.

During random walks, swimmer may exit the biofilm domain. When the swimmer reaches the domain boundary, a new swimmer is introduced with a velocity oriented towards the interior of the domain while the original trajectory is stopped at the boundary.

A local sensitivity analysis (Figure 1) is performed by comparing basal simulation obtained with $\gamma =\beta =ϵ=1$ (v₀ and v₁ where taken as in Appendix 1—table 3) with 3 simulations where $\gamma$, $\beta$ and $ϵ$ are alternatively set to 0, resulting in 3 alternative models where the speed or the direction selection or the random term is turned off. The interaction between the speed selection term (set by $\gamma$) and the random term is illustrated in Appendix 2—figure 3 where 5 repetitions of the same trajectory of a simplified Langevin equation (11) are displayed with or without friction ($\gamma =1$ or $\gamma =0$), but with the same random seed for the stochastic term so that the stochastic part is strictly identical.

To analyze the impacts of the non-dimensionalized swimming parameters $\gamma$, v₀, v₁, $\beta$, $ϵ$ on the locomotion behaviour, a global sensitivity analysis has been performed. The parameter space ${[0,1]}^5$ was uniformly sampled with n = 1000 points using the Fourier Amplitude Sensitivity Test (FAST) sampler of the SALib library that is the function SALib.sample.fast_sampler.sample (Cukier et al., 1973; Saltelli et al., 1999). We note that the interval $[0,1]$ covers a large parameter domain for some parameters, in particular $\beta$ which remains small after inference. For this parameter, the sensitivity analysis will show potential impact on the output, that may be ineffective in the parameter range of the inferred model.

For each point in the parameter space, a forward simulation is conducted on a population of swimmers on a representative biofilm extracted from the dataset (first batch of the B. pumilus dataset). Trajectory descriptors are then extracted and taken as observable of the sensitivity anaylsis that requires both the parameters sampling and the associated descriptors. Sobol indices of first order are then returned and pairwise partial correlations matrix has been calculated. Convergence of the Sobol indices has been checked by taking sub-samples containing less than $1,000$ points.

To allow inter-species comparisons in plots, the data and model outputs are re-normalized with common reference values $A_{ref}$ and $V_{ref}$ defined as the average of the species reference values (see Table 2 for values). Uni-dimensional distributions (Figure 5 upper panel, Figure 6b upper panel, Figure 7a, upper panel, and Figure 7b) were obtained with the gaussian_kde function of scipy.stats. T tests for mean comparison were performed using scipy.stats ttest_ind.

Two-dimensional distribution plots (Figures 5 and 6 b, Figure 7a lower panels) were obtained by first plotting the two-dimensional point cloud and approximating the point distribution with a gaussian KDE using scipy.stats gaussian_kde function. Then, the gaussian kde is evaluated at each point of the point cloud and quantiles 0.05, 0.5, and 0.95 of the resulting values are computed. Finally, quantile isovalues are plotted and the point cloud and the KDE are removed (see Appendix 4—figure 1 and Sec. KDE computation for details): this procedure ensures to enclose 5, 50% and 95% of the original points, centered in the densest zones of the initial point cloud.

Ternary plots (Appendix 2—figure 6) were obtained by first computing the contribution of each term of equation (4) to acceleration estimate. Namely, note

We compute the proportions $A{(k)}_i^s$ for $k\in \{b,\nabla b,\eta \}$,

Points $(A{(b)}_i^s,A{(\nabla b)}_i^s,A{(\eta )}_i^s)$ are then plotted in ternary plots using the Ternary python package (Weinstein et al., 2019) and approximated by gaussian KDE. Isolines are finally plotted as previously described.

To construct the plot in Appendix 1—figure 2, pairwise correlation of the biofilm density in the 384 samples has been computed (scikit-learn pairwise_distances, ‘correlation’ metric parameter Pedregosa et al., 2011), and the resulting similarity matrix has been displayed using Seaborn package clustermap function (Waskom, 2021) after hierarchical clustering (scipy.cluster.hierarchy linkage function Virtanen et al., 2020). Additional permanova has been computed to assess the significance of between-group dissimilarities using stats.distance package permanova function (scikit-bio development team, 2020 ).

All the image pre- and post-processing, calculations and statistics have been performed with custom scripts using the standard python libraries numpy (Harris et al., 2020), scipy (Virtanen et al., 2020), imageio (Klein, 2021), and pandas (McKinney, 2010). The forward swimming problem computation is computed using customed scripts built upon numpy (Harris et al., 2020) and H5py (https://www.h5py.org). Sensitivity analysis has been conducted with the SALib library (Cukier et al., 1973; Saltelli et al., 1999) (Sobol index, function SALib.analyze.fast.analyze) and the pingouin library (Vallat, 2018) (PCC, pcorr method). The Bayesian inference has been conducted using the STAN library (Stan Development Team, 2018) through its python interface pystan (Riddell et al., 2021). All plots have been made with the matplotlib python library (Hunter, 2007).

The whole python code have been made available and accessible at the following git repository https://forgemia.inra.fr/bioswimmers/swim-infer.

To illustrate the impact of each parameter on the interplay between the host biofilm and the swimmers trajectories, the model (2) was first computed on two mock biofilms. The first one is a square linear density gradient and the second is composed of large pores on a textured background mimicking the dense biofilm zones (Appendix 2—figure 1a). A basal simulation is computed with $\gamma =\beta =ϵ=1$ and will be used later on as reference for comparisons. These three parameters are alternatively set to zero to assess the resulting trajectories when the speed selection, the direction selection or the random term is shut down. Suppressing speed selection results in rectilinear trajectories ($\gamma =0$, Appendix 2—figure 1c), which is rather counter-intuitive since the remaining terms are designed to tune the direction. A discussion of this phenomenon is provided in Appendix 2 Influence of inference and stochastic terms on the trajectory descriptors. When suppressing direction selection ($\beta =0$, Appendix 2—figure 1d), the trajectories are no longer drifted downwards the gradient in the upper panel as in the basal simulation, and no longer follow the pores (lower panel). If the stochastic term is shut down ($ϵ=0$, Appendix 2—figure 1e), the trajectories directly go down the gradients and are trapped in the center of the image in the upper panel. When a pore is found along the run, the swimmer keeps following it without being able to escape the pore any longer unlike the basal situation (lower panel).

The link between the model parameters and the global trajectory descriptors introduced in Section Characterizing bacterial swimming in a biofilm matrix through image descriptors is less intuitive. A global sensitivity analysis of the trajectory descriptors (mean acceleration and speed, distance, displacement and visited areas) with respect to the parameters $\gamma$, v₀, v₁, $\beta$ and $ϵ$ is conducted in Model sensitivity analysis by computing their first order Sobol index (SI) and their pairwise correlation coefficient (PCC). The sensitivity analysis shows that the mean speed is mainly influenced by $\gamma$ and $ϵ$ with slightly negative and positive impact respectively, while acceleration is rather influenced by $\beta$ and $ϵ$ with positive impact. The link between the parameters and the other descriptors is more complex, including non linear effects (strong SI and small PCC) and parameter interactions (higher SI residuals, see Appendix 2—figure 2).

Appendix 2—figure 1

Download asset Open asset

Model sensitivity analysis

The link between the model parameters and the global trajectory descriptors introduced in Section Characterizing bacterial swimming in a biofilm matrix through image descriptors is not intuitive. A global sensitivity analysis of the trajectory descriptors (mean acceleration and speed, distance, displacement and visited areas) with respect to the parameters $\gamma$, v₀, v₁, $\beta$ and $ϵ$ is conducted by computing their first order Sobol index (SI) and their pairwise correlation coefficient (PCC).

Appendix 2—figure 2

Download asset Open asset

The residual variance is small for the median speed and acceleration but slightly larger for the distance, displacement and visited area indicating larger effects of parameter interactions for these outputs, that is output variations induced by joint shifts of the parameters (Appendix 2—figure 2). The SI of the parameters v₀ and v₁ are negligible, except for the displacement and the visited area. The parameters $\gamma$, $\beta$ and $ϵ$, that is the three weights associated to each component of the state Equation 2, are more influential. Distance and speed have several main drivers. The distance is impacted nearly equally by $\gamma$, $\beta$ and $ϵ$ and the PCC of these parameters is quite small, indicating that these parameters may induce indistinctly negative or positive variations of the travelled distance, except for $ϵ$ which is slightly negatively correlated. The median speed is mainly impacted by $ϵ$ (slightly positively) and $\gamma$ (slightly negatively), with relatively small PCC (Appendix 2—figure 2). The mean acceleration, the displacement and the visited area are preponderantly impacted by a main driver: the mean acceleration and the visited area are particularly impacted by $ϵ$, the stochastic term weight, with positive influence. The displacement is mainly influenced by $\gamma$ with no preponderant variation direction (null PCC, Appendix 2—figure 2).

Friction and random term in Langevin equations

To illustrate the interplay between the friction and the random term during a random walks, we solve the problem

(11) $\begin{array}{cc}\hfill \mathrm{d}v=-\underset{friction}{\underbrace{\gamma v\mathrm{dt}}}+\underset{\text{random term}}{\underbrace{\eta \mathrm{dt}}}& \end{array}$

(12) $\begin{array}{cc}\hfill v(0)=(0,0)& \end{array}$

(13) $\begin{array}{cc}\hfill X(0)=(0,0)& \end{array}$

in an unconstrained domain, with $\eta$ a 2 dimensional white noise with unitary variance. The friction parameter $\gamma$ is alternatively set to 1 (Appendix 2—figure 3, upper panel) or 0 (Appendix 2—figure 3, lower panel). We note that the random seed is the same for the simulations with or without the friction term, so that the stochastic contribution is completely identical in the upper and lower panels. The trajectories produced without the friction term are much more regular and rectilinear that those produced with the friction term, that are much chaotic.

The reason of that behaviour may come from the null mean of the white noise. Roughly speaking, in average, the acceleration shows small variations around zero which leads after temporal integration to regular speeds and rectilinear-like trajectories. By contrast, the friction term reduces the particle inertia, enhancing the impact of the stochastic term, which produces much more chaotic trajectories.

Appendix 2—figure 3

Download asset Open asset

Impact of the stochastic term

We illustrate the impact of the random walk term on the overall swimmer trajectory with Appendix 2—figure 4. In this figure, we display two trajectories computed from model (2) with identical parameters ($\alpha$, $\beta$, v₀, v₁, $\gamma$ and $ϵ$), initial condition, host biofilm and time length. Different random samplings of the stochastic term of Equation (2) lead to these very different trajectories. This example illustrate the difference between identifying population-wide characteristics and inferring true trajectories: while the later try to detect the differences between the two trajectories (i.e. in this example, identifying and smoothing the different stochastic samples leading to these trajectories), the former focuses on the common features between these apparently different trajectories.

Appendix 2—figure 4

Download asset Open asset

Influence of inference and stochastic terms on the trajectory descriptors

We wonder if the uncertainty sources involved in the inference process and in the stochastic term of the random walk have a decisive impact on the trajectory descriptors. To address this question, a first dataset is assembled by integrating in time Equation 2 for given parameters (see Appendix 1—table 3), initial conditions and host biofilm. Then, this dataset is used as inputs of the inference method to infer the initial parameters (ground truth). Another dataset is produced by replacing the initial parameters by the inferred parameters. We note that we take the same seed for the random number generator than for the initial dataset, so that the only uncertainty that has been introduced until this step comes from the inference procedure. Finally, we produce a last dataset by solving the model with the same inferred parameters as in the second dataset, but changing the seed of the random number generator. Hence, this last dataset involves uncertainties coming from the stochastic terms and from the inference process. This variation results in modifying the sampling of the stochastic terms and leads to strong modifications of the trajectories, like in Appendix 2—figure 4.

At end, the trajectory descriptors are computed and plotted in Appendix 2—figure 5. We can see that the trajectory descriptor distributions are very similar across the different dataset, except for the total distance and the displacement where discrepancies can be noted. However, these differences are relatively small compared to the mean and the width of the distributions. We can also observe that the interactions with the underlying biofilm is very well conserved, even when the sampling of the stochastic term is very different. This observation grounds the initial guess that these trajectory descriptors captures common global features of the different trajectories rather than specificities of given trajectories.

Appendix 2—figure 5

Download asset Open asset

Relative impact of the different swimming processes on the species swim

The ternary plot presented in Appendix 2—figure 6 shows the balance between the different swimming processes. The contribution of each term of equation (4) to acceleration estimate was first computed. Namely, the relative value of the speed selection term $\parallel s{(b)}_i^s\parallel$, the direction selection term $\parallel s{(\nabla b)}_i^s\parallel$ and the stochastic term $\parallel s{(\eta )}_i^s\parallel$ where

$s{(b)}_i^s=\parallel \gamma (v_0^s+b(t,X_i^s(t))(v_1^s-v_0^s)-\parallel V_i^s(t)\parallel )\frac{V_i^s(t)}{\parallel V_i^s(t)\parallel }\parallel ,$

$s{(\nabla b)}_i^s=\parallel {\beta }^s\frac{\nabla b(t,X_i^s(t))}{\parallel \nabla b(t,X_i^s(t))\parallel }\parallel ,\text{and}\mathit{\hspace{1em}}s{(\eta )}_i^s=\parallel {\eta }^s\parallel .$

The proportions $A{(k)}_i^s$ of each process was computed with

$A{(k)}_i^s={\frac{s{(k)}_i^s}{s{(b)}_i^s+s{(\nabla b)}_i^s+s(\eta )}}_i^s.$

for $k\in \{b,\nabla b,\eta \}$. As the contribution of the direction selection was limited compared to the other processes, we zoomed in the plot near the edge $\parallel s{(\nabla b)}_i^s\parallel =0$ to allow inter-species comparisons.

Appendix 2—figure 6

Download asset Open asset

Different inference models were designed and tested from the dimensionless state equation (2).

We illustrate the process of visualization of multiple point distributions in the same graph using KDE and isolines enclosing specific proportions of the data in Appendix 4—figure 1. A point cloud is first approximated with a Gaussian KDE. Then, the value of the gaussian KDE is evaluated in each point of the original point cloud, which allows to map the 2D map into a 1D set where order relation can be defined. Specific quantiles of the resulting values are computed (namely quantile 0.05, 0.5 and 0.95). By definition, the quantile 0.05 separate 5% of the points of the original dataset (the 5% lowest Gaussian KDE values) from the remainder of the data set. The isoline corresponding to the quantile 0.05 then also separates in the 2D map the 5% lowest Gaussian KDE values from the others.

Appendix 4—figure 1

Download asset Open asset