Bayesian inference of wall torques for active Brownian particles
Sascha Lambert Merle Duchene Stefan Klumpp


Interactive Surgical Liver Phantom for Cholecystectomy Training
    [
    September 9, 2024
===============================================================

















    




§ INTRODUCTION




Self-propelled particles are a widely studied type of active matter that includes a broad range of systems, from microbial motility to animal flocking <cit.>. Active Brownian Particles (ABP) form one of the simplest models of such systems, based on a fixed active propulsion force and diffusion to describe the motion. These minimal descriptions are often extended to include interactions between the particles <cit.>, with external fields<cit.>, or with their environment <cit.>. The latter is important as self-propelled particles often move in complex environments characterized by confining walls or interactions with obstacles. Many studies have investigated the microscopic origin of alignment with obstacles based on steric and/or hydrodynamic interactions<cit.>. Within a phenomenological description, such alignment can be implemented within ABP models as a 'wall torque', an empirical torque that provides the observed alignment with surfaces <cit.>. Experiments can inform some of the underlying parameters of such extended models, while others are strictly of an empirical nature, combining multiple physical phenomena.

Our goal in this study is two-fold. First, we demonstrate that the aforementioned empirical torque model is a quantitatively plausible description of an active rod's steric interaction with a surface and explore how its parametrization influences accuracy. For this, we utilize a classical least-squares approach to fitting the torque model to simulations of an active rod and evaluate the relevant residuals that describe the errors in the predicted dynamics. Second, we describe a procedure to learn the functional form of wall interactions, such as the wall torque from experimental data. This is done using a Bayesian inference approach, where the posterior model parameters are evaluated from trajectories of self-propelled particles. Notably, this approach does not require any information about the particle's orientation, which might be experimentally inaccessible.



§ MODEL


Our base model is an Active Brownian Particle (ABP)<cit.>  in two dimensions with anisotropic diffusion, obeying the following set of Langevin equations:

    ṙ(t)    = v_0 φ̂e_x + 1/k_BTφ̂^-1 D_Tφ̂F_wall + √(2)φ̂^-1√( D_T)ξ(t), 
    φ̇ (t)    = D_R/k_BTτ_ABP + √(2D_R)χ(t).

They describe the time evolution of the active particle's location 𝐫 and its orientation φ using the active velocity v_0 and the rotational diffusion D_R. In contrast to the simplest ABP models, we allow for anisotropic translational diffusion, as we will use this model to describe elongated (rod-shaped) particles explicitly. The translational diffusion is thus described by a diagonal diffusion matrix D_T=diag(D_T^∥, D_T^) in the particle's coordinate system. We use φ̂ to denote the rotation matrix linking the lab's coordinate system to the particle's coordinate system. The diffusion is set to values representing a spherocylindrical rod, using the parametrization of 
Lüders et. al. <cit.>:


    D_T^∥   = k_BT/2 πη L(log (p) - 0.1404 + 1.034/p - 0.228/p^2), 
    
            D_T^   = k_BT/4 πη L(log (p) +0.8369 + 0.5551/p - 0.06066/p^2), 
    
            D_R    = 3 k_BT/πη L^3(log (p) - 0.3512 + 0.7804/p - 0.09801/p^2).

p is the aspect ratio of the spherocylinder, and L is its length. This parametrization is valid for p=1-30 <cit.>.


The ABP may interact with walls, which exert force and torque on the particle.
The term F_wall=-∇ V_WCA in eq. <ref> is the volume exclusion force exerted by walls. We use WCA repulsion <cit.>, a Lennard-Jones potential with a truncated attractive region:

    V_WCA(x) =
         
        V_LJ(x) - V_LJ(x_c)    for  x < x_c =2^1/6/2σ, 
    
        0    for  x ≥ x_c,

with

    V_LJ(x) = 4ϵ[ ( σ/x)^12 - ( σ/x)^6].

Here, σ is the length scale of the interaction. For our simulations, we set the hardness to ϵ=4k_BT.

In addition, the ABP gets subjected to an empirical torque τ_ABP, factorized into a force term proportional to the instantaneous repulsion exerted by the wall (which depends on the distance from the wall) and an angle term dependent on the ABP's incidence angle with the wall:

    τ_ABP = |F_wall|· f_ABP(θ_wall - φ).

The angle term f is an arbitrary function of the angle of incidence between the swimming direction φ and the wall's normal direction θ_wall. We aim to learn this function using information gathered from observations of an active rod's trajectory. For this, we parametrize the unknown function in terms of its spectral components

    f_ABP(φ) = ∑_n=1^N_Pα_n
        ·sin(n[ θ_wall - φ]).

up to order N_P. All anti-symmetric cosine terms are suppressed by demanding symmetry with respect to the wall-normal, leaving only the sine terms. The factors α_i determine the modes' strengths. The first mode (n=1 with α_1<0) always produces a torque away from the wall, irrespective of the angle of incidence. This creates an additional repulsion mechanism that leads to trajectories moving away from the wall in addition to the WCA repulsion. The second mode (n=2) creates a torque that aligns the ABP with the wall tangent. When the particle is heading toward the wall, the torque turns it away from the wall, whereas when the particle is headed away from the wall, the torque turns it back to the wall tangent. With additional higher frequencies, more complex torque functions can be described.



 §.§ Training Data generation


The Bayesian method outlined in this paper provides a framework for learning the wall torque function f from experimental data on a swimmer's location over time. However, for the context of this paper, we only employ synthetic data, as this allows for better diagnostics of the method. The data is generated using model variants described in the previous section. We generate two types of trajectories:

    
  * Type-A: ABP trajectories with predefined wall torque functions f. These provide ground truth data to validate the Bayesian inference method.
    
  * Type-R: Rod trajectories from a variant of our ABP model describing rod shape particles explicitly. In this case, the empirical wall torque is disabled, f=0. To capture steric interactions of a rod, the ABP gets equipped with N_S test sites equidistantly distributed at locations ± piσ/2N_S in the front and back of the ABP. The WCA repulsion forces F_± sat, i are evaluated at the locations of the test sites and produce torques piσ/2N_S· F_± sat, i×φ̂ e_x.

Both types of simulations are set up like the example in fig. <ref>. The ABP and the rods are initialized outside the WCA interaction region in front of the wall. The initial orientations were chosen on a grid from -60^∘ to 60^∘. Integration is performed using an Euler-Maruyama discretization and terminated when the particle is sufficiently separated from the wall. 





§ LEAST-SQUARES ANALYSIS


To evaluate the performance of the empirical torque model in describing the steric interactions of rods, we've extracted the dynamic torque from type-R simulations.

We generated 2500 rod trajectories (100 per initial orientation), each for a range of parameter sets. The rate of change φ̇(t) ≈ (φ(t+h) - φ(t))/h in the orientations, conditionally averaged for given distance and orientation, can then be used to calculate the effective torque on the rod:

    τ(t) = k_BT/D_rot ⟨φ̇(t)⟩.

Fig. <ref>a shows an example of the resulting torque data, parametrized by the distance to the wall and the rod orientation. The system's dynamics primarily occur in a crescent moon-shaped region, where the rod is inbound to the wall. The white regions are not sampled in the simulations, i.e., they are typically not visited during the scattering process. When the rod is pointed away from the wall, the torque quickly vanishes as the active motion and the repulsive interaction move it out of the interaction zone of the WCA potential.



We then fit eq. <ref> to the data by minimizing the squared residuals:

    (α_n)_opt =(α_n)min∑_t,k[ τ^(α_n)(φ^k(t), d^k(t)) - τ^k(t)]^2,

where the index k runs over the simulated trajectories.
Assuming Gaussian errors in the calculated torques (which is expected given the diffusive structure of the equations of motion), this gives us a point estimate of the most likely torque parametrization. We calculate error bands using 68% confidence intervals obtained from bootstrapping <cit.>. Fig. <ref>B shows the best fit obtained for the leading mode (n=1). The corresponding amplitude is  (α_1)_opt=2.03σ±0.08σ.

We've tested the influence of several parameters on the model. Variations in the number of test sites for the repulsive interaction with the wall, the Péclet number of the active rod, and the wall hardness did not significantly influence the performance of the empirical model or the parametrization of f (see Supplementary Material).





In contrast, the rod aspect ratio p, the mode frequency n, and the mode count N_P affect the fit significantly. Fig. <ref>A and B show the fit results and performance for a single-mode description of f. Note that we've varied n continuously instead of only considering integer values. We find a linear relationship between 1/n and its optimal amplitude (α_n)_opt. This suggests that the rod dynamics at small angles of incidence, where f∝α_nn, are of primary importance to the structure of f.
Furthermore, we see a linear relationship between the amplitude and the rod aspect ratio, corresponding to a modulation of the lever arm length. Surprisingly, the tangentially aligning mode (n=2) is not the optimal configuration. Instead, we found that a slightly bigger n≃ 1.1 shows lower residual errors for all aspect ratios, as shown in fig. <ref>C.

When using multiple modes (discrete n=1,2,...), we find very similar results to the single-mode case in terms of performance and dependency on the aspect ratio. The ratios of optimal mode amplitudes show a non-linear dependency on p (see Supplementary Material).  
Using up to four modes improves the model's performance, albeit by less than 1% compared to the single-mode model. We did not find that using more than four modes leads to any further improvement of the model (see fig. <ref>D).

From these results, we conclude that using the empirical torque function f successfully introduces a mechanism that captures the steric repulsion seen in rod simulations. However, different modes of the torque function result in approximately equal performance, indicating that the choice of torque function is not unique and that modes are somewhat interchangeable. Likewise, extending the representation to include higher-frequency components does not significantly improve the model quality after four modes. 



§ BAYESIAN ANALYSIS


We now turn to a Bayesian approach <cit.> to learn the torque function f from observations. This approach provides intuitive and easily interpretable <cit.> results that quantify our knowledge about f that is gained by observing the system dynamics. Specifically, we construct the posterior probability distribution of f. 



 §.§ Method Outline


Bayesian inference targets the posterior distribution p(θ|Y), where θ denotes the set of parameters of interest, in our case, the parametrization (α_n) of the torque function f, and Y denotes the data used for the inference, here consisting of sequences of x- and y-positions that form the trajectories of the rod (but not its orientation). We use the tip of the rod as the measured location, as it is the primary interaction point of the rod with the wall. The posterior distribution is calculated from Bayes' rule p(θ|Y) ∝ p(Y|θ)p(θ), using the likelihood p(Y|θ) and the prior p(θ). We use an uninformative uniform prior for all components of θ, only limiting the interval so that the ABP does not turn more than 90^∘ in a single integration time step. Furthermore, we assume a flat prior for all unobserved orientational states.

The likelihood p(Y|θ) is computed for the (type-A) ABP model. Here, we exploit that our observation data originates from another simulation (type-A for test cases, then eventually type-R) with a finite integration time step. We adopt the same time steps for the ABP, thus avoiding introducing another level of discretization error. If the method is applied to real experimental data, one will introduce a discretization error corresponding to the data frequency (for example, the framerate of the experimental video). Using a single time step between observations conveniently decouples the rotational and translational diffusion of the ABP, making it possible to calculate the likelihood just from the particle's location while marginalizing the ABP's orientation with a filtering approach.

After discretization with an Euler-Maruyama scheme, equation <ref> becomes

    r_i+1 - r_i - v_0 φ̂_i e_xΔ t - 1/k_BTφ̂^-1_i D̂_Tφ̂_i F_wallΔ t
    
            =√(2)φ̂_i √(D̂_TΔ t)w_i(t) ≡o_i^φ_i,

yielding the displacement the particle experiences through translational diffusion as the offset o_i^φ_i. The random numbers w_i are normal-distributed. Therefore, the offset is normal-distributed as well:

    o_i^φ_i∼𝒩(0,2φ̂_iD̂_Tφ̂_i^-1Δ t).

Then, the joint likelihood of all the offsets, conditioned on the orientations φ_i of the ABP at each time step, is

    p(o_1^φ_1...o_N^φ_N|θ, φ_1...φ_N)=∏_i^N𝒩(o_i^φ_i|0,2φ̂_iD̂_Tφ̂_i^-1Δ t)




 §.§ State Marginalization


This equation needs to be marginalized for all orientations φ_i, which is not analytically possible due to the trigonometric functions in φ̂_i. As we're dealing with a state space model, an effective strategy is using a sequential filter, which calculates the posterior of φ_i given all previous states φ_1...φ_i-1. The filter then steps through all time steps in sequence to construct the joint probability distribution p(φ_1...φ_N|θ) iteratively. We considered a range of commonly used filters and found that particle filters <cit.>  work very well, unlike linear Kalman filters. We also decided against non-linear extensions of the Kalman filter as the particle filter produces an unbiased estimation of the posterior density, which allows us to sample from the true posterior via pseudo-Marginal Metropolis-Hastings (PMMH) methods <cit.>.

The particle filter approximates the posterior distribution of the rotational state φ_i as a sum over k=1...N_p delta distributions at locations φ̃_i^k, the name-giving 'particles':

    p(φ_i | θ, φ_1...φ_i-1) ≈1/N_p∑_k=1^N_pδ(φ_i - φ̃^k_i-1)

Filter particles, the discretization samples for the orientational states,  are characterized only by the unobserved orientational degree of freedom of ABPs, but, together with the positional data, can be used to sample full ABP trajectories that include the unobserved latent orientations. They follow the orientation dynamics of eq. <ref> with the torques that follow from the position data.
With this, we can marginalize equation <ref> as

    p(Y|θ)    ∝∫ p(o_1^φ_1...o_N^φ_N|θ, φ_1...φ_N)p(φ_1...φ_N)φ_1...φ_N 
       ≈∏_i=0^N1/N_p∑_k=1^N_p p(o_i^φ̃_i^k|θ, φ_1^k...φ_i-1^k)

Note that the offsets are the only stochastic terms in Y when conditioned on the latent orientations, which is why, inside the integral, the offsets form the entirety of the observed data. The particles' initial states are sampled from the prior φ̃^k_0 ∼ p(φ_1). After that, the particles follow the orientational dynamics of equation <ref> and are resampled with weights w_i^k=p(o_i^φ̃_i^k|θ, φ_1^k...φ_i-1^k)p(φ_i-1) after each time step. We do the resampling step before applying rotational diffusion to reduce particle degeneracy, where the orientation distribution is poorly represented as many identically resampled particles.

To summarize, we employ the following procedure to construct an estimation of the posterior density of the parameter set θ:


    
  * Initialize N_p filter particles (discrete orientation values) sampled from the prior p(φ_0).
    
  * Likelihood estimation: Calculate the offsets o_i^φ̃_i^k for all particles and update the partial product in equation <ref>.
    
  * Assimilation: Resample the particles based on the latest likelihood estimation and the priors on the orientations.
    
  * Dynamics: Apply the equation of motion <ref> to all filter particles, i.e. update the orientation by rotational diffusion and the torques.
    
  * Repeat from 2 until all data has been used.




 §.§ Posterior Density calculation


The procedure described above does not normalize the posterior density, as the marginal probability is intractable to calculate. We address this with two different techniques. The first is to use a Markov-Chain Monte-Carlo (MCMC) method, the PMHH algorithm <cit.>,  to sample directly from the posterior. This method relies on the particle filter as an unbiased estimator of the posterior density. For visualizations of the posterior, we used the second approach, discretizing the parameter space θ and doing an exhaustive sweep. The discretization is then readily normalized.





 §.§ Posterior Envelope

To tune the inference performance in large parameter spaces, we employ a procedure to adapt a multivariate Gaussian envelope to the posterior. This can be done with relatively few filter evaluations, and the resulting distribution is then used directly as a proposal for the Markov chain.

Our fitting procedure produces a sequence of Gaussians 𝒩(θ_j, Σ_j). Starting with an initial guess θ_0, Σ_0, we draw k=1...N_T parameter vectors from a test distribution θ_j^k ∼𝒩(θ_j, zΣ_j), with the windows size z serving as a scaling of the search window. We then estimate the posterior densities p(Y|θ_j^k) at the test points using the particle filter and calculate the next iteration as

    θ_j+1   = Σ_k=1^N_Tw_j^kθ_j^k
    Σ_j+1   = ∑_k=1^N_T w_j^k (θ_j^k - θ_j+1)^T (θ_j^k - θ_j+1),

with weights

    w̃_j^k = p(Y|θ_j^k)/𝒩(θ_j^k|θ_j, zΣ_j)   w_j^k = w̃_j^k/∑_i=1^N_Tw̃_j^k.


This weighs test samples based on the targeted posterior density and compensates for the test distribution, which biases the sampling towards its mean θ_j. This procedure is iterated until convergence is observed in the covariance matrix. Convergence can be monitored from the entries of Σ_j (see Supplementary Material). We've found that N_T=1024 test samples, 1500 particles, and a search radius of z=1.5 work well in our scenarios, but we have not tuned the procedure in detail.





 §.§ B-Spline representation


In our tests, we've found that using the spectral representation from eq. <ref> is causing numerical issues when using more than 3-4 modes. In particular, the posterior mass collapses to a very thin distribution that is not aligned with any specific α_i dimension. The strong correlation makes it numerically difficult to find and represent the envelope of the posterior accurately, given the noisy density estimates from the filter. We attribute this to the fact that the Fourier modes are highly non-local and are affected by all observations.  By switching to a quadratic B-spline representation, we can alleviate the problem. The B-spline is of the form

    B^(2)(x)    = ∑_i=0^N_Bα_i·κ((N-3)x + i - 1/2)
     with κ(u)   = 1/2u^2                             0≤ u ≤ 1,
    1/2(-2u^2 + 6u - 3)    1≤ u ≤ 2,
    1/2(3 - u)^2      2≤ u ≤ 3, 
    
    0    else

and represents functions in the base system κ(.), stitched together from quadratic functions. The most important property of this spline for our use case is that κ has compact support and thus is not affected by all observations, like the spectral representation.




 §.§ Validation on type-A data

We tested the procedure on data generated from (Type-A) ABP simulations with known torque function f. An example is shown in fig. <ref>A and B. We used a 2-mode torque (with ground truth α_1 = α_2 = 10σ) and generated 20 trajectories. We then inferred the amplitudes to check if the implementation correctly recovers the ground truth. The posterior density shows a highly complex structure that depends on the particular observations used for the inference. Most of the probability mass concentrates into a Gaussian near the ground truth. In all our tests, the 99% highest density interval (HDI) reliably captures the ground truth. One can see a clear anticorrelation between the amplitudes of the two modes, indicating that the modes are, to some extent, interchangeable.

We also tested the inference of a B-spline-represented torque function on the same data (fig. <ref>C). This representation also reliably captures the torque function, at least for inbound trajectories. Strikingly, the torque function is not well-learned in the outbound orientations when the particle moves away from the wall. This matches the observations from the least-squares approach, where dynamical data is collected in a crescent-shaped region of the φ-d-space. Information about the torque function for a specific orientation is primarily gained when the particle is inside the WCA interaction region and pointing in the respective direction. As the center of the crescent shape is never realized in the dynamics, no information about the torque can be gained.



 §.§ Inference on rod data




Finally, we learn the torque function f from 20 Type-R simulations of active rods, which serve as a proxy for experimentally collected data. As these simulations now use true steric interactions with the wall instead of a predefined torque function, no ground truth for f is known. Results for rods with aspect ratio p=1.5 are shown in fig. <ref>. Notably, the posterior density coincides with our fitting results for trajectories that are inbound to the wall, giving credence to the previous interpretation that this is the regime where the torque function is most important. The Bayesian inference furthermore reveals that most knowledge about f is gathered in the inbound direction. The HDI remains broad for parallel and outbound trajectories.

Model selection imposes strong constraints on the torque function: Using two modes limits the bandwidth of the function, correlating the inbound and the outbound direction (fig. <ref>A). In this case, the HDI remains finite and relatively small (compared to the prior) for outbound movement, as knowledge about the regime is learned from information that is gathered from inbound movement.

Using a higher resolution representation of f, such as the B-splines representation shown in fig. <ref>B lifts these correlations. Then, the posterior f densities for inbound and outbound movement are constructed from different observations. This reveals two important properties. The first is that f is not strongly defined for small angles of incidence. While the inference strongly prefers negative torques (which corresponds to steric repulsion), there's no dominant magnitude that is needed to explain the trajectories. This may be explained by rotational diffusion being more important when a normally oriented rod experiences little torque, in combination with a lack of data on such rods, as they are in an unstable configuration. Secondly, the uncorrelated B-spline representation shows that the torque function does not impact the trajectory prediction in the outbound regime. The posterior coincides with the prior in this regime.







§ CONCLUSIONS


In this study, we have validated an empirical torque model to describe the steric interactions of active rods with walls. To this extent, we've simulated spherocylindrical active rods scattering with a wall and compared the resulting torques to an Active Brownian Particle model, which experiences a torque described by an orientation-dependent function f. By employing a least-squares approach, we demonstrated how the ABP model accurately fits the simulated rod trajectories. We've looked at the functional structure of f and found that it is primarily influenced by the rod's aspect ratio. We've discovered no significant impact of the discretization of the rod or its activity on the model's accuracy.

We also introduced a Bayesian inference framework to learn the torque model from observational data. This method works robustly in our test cases, where we generated trajectories with predefined torque functions, which the inference could reasonably recover. The posterior density of the torque function directly highlights where information about the system is gained and where the model stays elusive. We've shown that care must be taken when selecting the representation of f, as insufficient complexity introduces unwarranted correlations that suppress details of the inference. We've applied the method to data obtained from rod simulations and found that the empirical torque model is effective for rods inbound to a wall, whereas the outbound direction is, in practice, not affected by the choice of torque model.

The Bayesian inference method presented here relies on the assumption that the system dynamics can be described by a single Euler-Maruyama step between observation times. While this is exact for our synthetic data, it will be an approximation when the method is applied to gather information from real experiments, and its accuracy remains to be explored. Finally, we want to emphasize that the Bayesian approach described here to infer wall torques is general and also applicable to inferring other properties of active particles. 

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project ID 446142122. 
Simulations were run on the GoeGrid cluster at the University of Göttingen, which is supported by DFG (project IDs 436382789; 493420525) and MWK Niedersachsen (grant no. 45-10-19-F-02).


eplbib


