Large-Scale Dynamics of Self-propelled Particles Moving Through Obstacles: Model Derivation and Pattern Formation

We model and study the patterns created through the interaction of collectively moving self-propelled particles (SPPs) and elastically tethered obstacles. Simulations of an individual-based model reveal at least three distinct large-scale patterns: travelling bands, trails and moving clusters. This motivates the derivation of a macroscopic partial differential equations model for the interactions between the self-propelled particles and the obstacles, for which we assume large tether stiffness. The result is a coupled system of nonlinear, non-local partial differential equations. Linear stability analysis shows that patterning is expected if the interactions are strong enough and allows for the predictions of pattern size from model parameters. The macroscopic equations reveal that the obstacle interactions induce short-ranged SPP aggregation, irrespective of whether obstacles and SPPs are attractive or repulsive. Electronic supplementary material The online version of this article (10.1007/s11538-020-00805-z) contains supplementary material, which is available to authorized users.

of collective dynamics studies what happens when a large number of agents, which can be animals, people, micro-organisms, crystals, etc., interact with each other. A particular focus is the emergence of large-scale order or patterns. Famous examples include global alignment in crystals (de Gennes and Prost 1993), lane formation for people (Feliciani and Nishinari 2016), waves and aggregation in bacteria (Shimkets 1990;Ben-Jacob et al. 2000), milling in schools of fish (Shaw 1978) or swarming in birds (Cavagna et al. 2010). All these examples have in common that local, small-scale interaction rules between individuals lead to global, large-scale patterns. These patterns are typically hard or impossible to predict from the local interaction rules; hence, their understanding requires the use of either extensive simulations or mathematical analysis.
Combining Collective Dynamics and Environmental Effects In many systems, one also needs to take into account the environment to be able to explain observed patterns in collective phenomena Jabbarzadeh et al. 2014;Park et al. 2008). For cells moving through a tissue, this environment often includes fibres and other components. For instance, it has been observed that many cell types have a tendency to move up stiffness gradients, a phenomenon termed durotaxis (Lo et al. 2000). In some of these instances, the effect on the substrate is negligible or the environment forms confining barriers affecting organism behaviour (Noselli et al. 2019). However, in many applications the interaction modifies the environment (either permanently or transiently) in a way that affects subsequent interactions. An example is the degradation of the extracellular matrix (ECM) caused by migrating cells (Baricos et al. 1995), which affects the ECM structure and hence future migration. In this work, we want to combine collectivity and environmental interactions and study the resulting patterns. Known examples of patterns created include travelling bands of large swarms of scavenging locust (Buhl et al. 2006;Topaz et al. 2008), the formation of paths in grass-land by active walkers (Helbing et al. 1997;Lam 1995) or aggregation of individuals (Bernoff and Topaz 2016). For metastasising cancer cells, it was observed that the invasion success depends on whether they move individually or as small clusters (Cheung and Ewald 2016).

Obstacles Can Emulate Complex Environments
In this work, we focus on a particular type of environment in which objects interact with moving particles. Example applications include pedestrians avoiding obstacles (Helbing et al. 2005) or animal herds or fish swarms moving through vegetation. Our focus, however, lies on micro-scale applications, such as pathogens moving through visco-elastic tissue (Celli et al. 2009;Harman et al. 2012) or immune cells migrating through fibrous ECM (Baricos et al. 1995). The importance of the environment is particularly true for sperm dynamics, where the surrounding fluid plays a key role in the emergence of collective motion. For example, clustering and large-scale swirling was observed in simulations of collectively moving sperm in Schoeller and Keaveny (2018) and Sokolov et al. (2007). In Degond et al. (2019), a model was proposed that couples the Vicsek model for collective dynamics with Stokes equations for a viscous fluid. However, sperm dynamics takes place in a complex fluid, whose constitutive properties cannot be characterised solely by a viscosity. In Tung et al. (2017), it was reported that sperm moving through a visco-elastic fluid forms small clusters, a behaviour not observed in a purely viscous environment. To approximate the complex environment, the introduction of immersed obstacles has been proposed (Kamal and Keaveny 2018;Majmudar et al. 2012;Wróbel et al. 2016). For example, Kamal and Keaveny (2018) the authors propose a model in which an undulatory swimmer swims in a fluid filled with elastically tethered obstacles, however effects of collective dynamics, i.e. multiple swimmers, were not investigated. In this article, we present a model for collective motion in an environment filled with spheres tethered to fixed points in space via linear springs, that play the role of obstacles. We will study the impact of this obstacle-based environment on the collective dynamics for a large number of self-propelled particles (SPPs).
Consequences for Aggregation Aggregation of populations of individual animals or cells has been observed in many contexts and is usually attributed to individuals' attraction (Parrish and Edelstein-Keshet 1999). A major finding of this work is that an initially homogeneous elastic environment can lead to particle aggregation in the absence of explicit attractive interactions between the particles themselves. Furthermore, we show that the obstacles can induce aggregation irrespective of whether they repel or attract the particles. This poses the question whether past conclusions about the cause of biological aggregation need to be reviewed: Some of the observed aggregative effect might have been caused by so far under-appreciated environmental interactions, rather than by direct interactions between individuals.
Individual versus Continuum Description From a modelling perspective, two approaches are common (Mogilner and Manhart 2016). One can formulate a system of individual-based models (IBMs), also called agent-based models, where the behaviour of each individual is assumed to be governed by separate, often stochastic ordinary differential equations. This approach has the advantage that the translation of modelling assumptions of the individual level is relatively straight-forward. However, few analytical tools are available to study IBMs and even if the system exhibits the desired property, limited insight can be gained as to why it does so. On the other hand, one can formulate a partial differential equation (PDE) model for the macroscopic quantities of interest, e.g. the space and time-dependent density of agents. A rich mathematical toolbox exists for the analysis of PDEs, which includes linear stability analysis, constructions of steady states as well as efficient simulation tools. Substantial progress has been made to establish systematic links between IBMs and the corresponding PDEs (Degond and Motsch 2008;Ha and Tadmor 2008). This allows to combine the advantages of both methods: straight forward translation of biological assumptions into the IBM and strong analytical tools for the PDE model. The self-organised hydrodynamics (SOH) approach (Degond and Motsch 2008) used in this work has been successfully applied, e.g. to fibre interactions (Peurichard 2016), bacterial swarms (Degond et al. 2018), sperm fertility (Creppy et al. 2016) or ant trail formation (Boissard et al. 2013).
Paper Structure In Sect. 2, we present the individual-based model, at whose basis lies the famous (Vicsek et al. 1995) model. This model describes SPPs that align their orientation with neighbouring particles, to which we add a short-ranged repulsion term. The environment consists of obstacles which are tethered via linear springs to anchor points fixed in space. SPPs and obstacles exert either repulsive or attractive forces on each other. Simulations of the IBM reveal the richness of possible patterns for this simple system, which includes clustering, trail formation and travelling bands, and

Formulation of the IBM
The following model describes self-propelled individuals interacting with obstacles. Biologically, the particles in our model typically represent cells (such as cancer, sperm or bacterial cells), but can also describe animals or people. We assume they align and repel each other. Alignment might be an active behaviour (e.g. birds might try and adjust their movement direction to each other), or a passive effect caused, e.g. by the individuals' shape (rod shaped bacteria might align upon collision (Peruani et al. 2006). Obstacles are elastically tethered to fixed anchor points. They mimic a complex environment, which acts on and reacts to particles, e.g. by repulsion. They can represent, e.g. a fibrous network which is relatively fixed in space, but whose components can be pushed upon interaction, after which they relax back towards their original position. Other applications could include plants such as trees or sea grass. In the following, we mathematically formalise these assumptions. The starting point for our investigation is an individual-based model (IBM), in which the dynamics of each component is described by individual equations coupled through interaction terms. We couple the famous Vicsek model for collective movement of self-propelled particles (SPPs) (Vicsek et al. 1995) with an environmental model, described by elastically tethered obstacles. Our IBM is set in n-dimensional space, where n = 1, 2 or 3. The two components and interactions are depicted schematically in Fig. 1. Several applications of collective movement, in particular when applied to cells, take place at the micro-scale. These regimes are typically friction dominated with negligible inertia (also called over-damped regime). We therefore formulate our model in this friction dominated regime.

Model Components
We model the following two types of agents: -Obstacles We consider a set of N mobile obstacles with positions X i (t) ∈ R n for i = 1, 2, . . . , N and time t ≥ 0. Each obstacle is tethered to a fixed anchor point Y i ∈ R n through a Hookean spring with stiffness constant κ > 0 and experiences friction with the environment with friction constant η > 0. -SPPs We denote by Z k (t) ∈ R n the positions of the k-the SPP at time t ≥ 0 for k = 1, 2, . . . , M. Each SPP has a body orientation α k (t) ∈ S n−1 and a selfpropulsion speed u 0 in direction α k . SPPs experience friction with the environment with friction constant ζ > 0.
Interactions We consider the following interactions: -SPP alignment We assume each SPP aligns its body orientation α k to the mean orientationᾱ k of body directions of SPPs in its neighbourhood with radius r A . This happens with an alignment frequency ν > 0 and is analogous to the famous Vicsek model for collective swarming (Vicsek et al. 1995). -SPP repulsion SPPs repel each other at short distances, which models sizeexclusion effects. Following Degond et al. (2015), we model this by an even pushing potential ψ : R n → R with typical spatial scale r R > 0. The force felt between two SPPs positioned at Z i and Z j is then given by ∇ψ Z i − Z j . -Obstacle-SPP interaction We assume the obstacles and SPPs exert a force on each other, which depends on the distance between them. Similar to the SPP repulsion, we describe this by an even interaction potential φ : R n → R with typical scale r I , yielding the force ∇φ (Z − X ) for a SPP at position Z and an obstacle at position X . In general, we assume this force to be repulsive; however, we will discuss the effect of an attractive force in Sect. 4.

Stochasticity
We include two sources of uncertainty, both modelled by independent Brownian motions: Stochastic effects in the obstacle position (with intensity d o ) as well stochastic effects in the SPP orientation (intensity d s ).

Model Equations
The effects described above can be modelled through the following coupled, stochastic ODEs. Note that in the absence of obstacles, the equations reduce to the time-continuous Vicsek model, described, e.g. in Vicsek et al. (1995). From here on, we work with the non-dimensional variables (but keeping the same names as introduced above), in particular we have chosen the domain size L as reference length and L/u 0 as reference time. The latter can be interpreted as the time it takes a freely moving SPP to cross the domain. We then obtain: where the mean directionᾱ k is defined via the mean flux J k bȳ The tether positions Y i are given and do not change in time. The operator P α ⊥ k in (1c) is an orthogonal projection onto α ⊥ k and ensures that if α k (0) ∈ S n−1 , then α k (t) ∈ S n−1 for all time. Note that we have scaled the interaction terms by the number of SPPs or obstacles to prepare for the kinetic limit of Sect. 3.1.

Remark 1 (Modelling choices)
In an attempt to create a minimal model, we did not include a number of effects. For example, as opposed to Degond et al. (2015), we don't model relaxation of the SPP orientation to the SPP velocity. Notice also that we did include repulsion between SPPs, but not repulsion between the obstacles. The former helps avoid collapse of the SPP density. For the obstacles on the other hand, this seems to be less likely due to their tethering in space. Also, there is no coupling to a surrounding fluid, which will be subject of future work.

Simulations of the 2D IBM
We simulate the IBM (1) in two space dimensions. In this work, instead of doing a more thorough investigation, we want to showcase what types of patterns can be created based on the environmental interactions, emphasising the need for a PDEbased description.  We distribute the fixed anchor points Y i using a uniform distribution on B and initialise the obstacle positions with X i (0) = Y i . Initial SPP positions Z k (0) and orientations α k (0) = (cos(ϕ k ), sin(ϕ k )) are both chosen at random with uniform distributions on B for Z k (0) and on [0, 2π ] for ϕ k . For the interaction potentials, we use kernels of the following shape where H (x) is the Heaviside function. These kernels are compactly supported on balls with radius r I and r R , respectively, and chosen to yield a continuous pushing force decreasing linearly. They are normalised such that the force mass is A I and A R , respectively. For simplicity, we choose all interaction radii to be the same, i.e. r I = r A = r R . We leave the following parameters constant: d o = 0, η = 1, d s = 0.1, A I = 1. We're left with five parameters: κ, ζ , ν, r I and A R . Figure 2 shows examples of the different patterns produced by different choices of parameters, and Fig. 3 shows some associated statistics. Corresponding videos can be found in Supp. Mat. The second row in Fig. 3 shows that in all three cases SPPs globally align, i.e. the variance in SPP direction decreases. We call the observed patterns: Moving clusters (κ = 100, ζ = 1, ν = 10, r I = 0.05, A R = 0.01), Trails (κ = 2.5, ζ = 10, ν = 100, r I = 0.15, A R = 0.002) and
Moving Clusters In this regime, tether stiffness and SPP-obstacle repulsion are high. The SPPs form very high density groups moving through the obstacles, whose displacement from the anchor points is relatively low. In Fig. 2a, we see how a larger cluster is split into two due to the obstacles, suggesting that the cluster size is controlled by the dynamics. This might also be the reason for the relatively large changes in mean SPP density over time seen in Fig. 3. Nevertheless, this pattern seems to be stable.
Trails Here, SPP alignment is strong, with low tether stiffness. The SPPs form stripes parallel to their movement direction, which at t = 10 seem to be very regularly spaced. Within the stripes, the SPPs are close together and consequently push the obstacles away from the trails, leading to large obstacle displacements. Interestingly, the trails become unstable and by t = 60, the SPPs form moving groups. We see this instability building up and the trails falling apart around t = 26 in Fig. 3b. The enlargements in Fig. 2 indicate that the instability of the trails might stem from the fact that the obstacles are not symmetrically displaced to the right and left of the moving trails.
Travelling Bands In this pattern, the spring strength is high and obstacle displacements are consequently small. SPPs now form bands normal to their direction of movement. At t = 60, we see in Fig. 2 that there appears to be a typical spacing between the bands. These patterns seem to be stable.

Obstacles Reinforce and Diversify Patterns
To assess the influence of the environment on the pattern formation, we compare to simulations of the model without obstacles, i.e. pure Vicsek-type dynamics with small SPP repulsion. In the inset in the second column in Fig. 2, we see that in all three examples there is no patterning in absence of the obstacles. Figure 3 shows that the alignment behaviour seems unaltered by the obstacles, however for moving clusters and trails the obstacles lead to much higher SPP densities. It is known that for some specific ranges of parameters, clusters and bands already appear in simulations of the Vicsek model alone (Vicsek et al. 1995). However, in the presence of obstacles their qualitative behaviour is different: the environment seems to reinforce such structures and the travelling bands appear to be regularly spaced, which is not the case for bands in the Vicsek model alone. In addition, the homogeneous phase (common to the Vicsek dynamics) appears to be less common here. Finally, we observe that also a completely new pattern emerges: trails.
The Need for a PDE Description The three patterns found by simulating the IBM show that the interactions between SPPs and obstacles can lead to a rich repertoire of patterns such as clustering, trails and travelling bands. While the system is relatively simple, the number of parameters (about 15) make it prohibitively expensive to explore fully the complete parameter space: Total computation time with an 1.80GHz Intel Core i7 CPU was on the order of magnitudes of hours for one simulation shown in Fig. 2. Simulation time increases with the number of SPPs and obstacles, with larger interaction radii as well as with larger forces produced by the dynamics (e.g., by  Fig. 2 for moving clusters (a), trails (b) and travelling bands (c). Solid and dashed lines mark averages for simulations with and without obstacles, respectively, shaded areas and dotted lines corresponding averages ± standard deviations. SPP densities are calculated for each SPP by calculating the density within a disc of radius r A = r R and dividing by the mean density in the domain clustering), which necessitates smaller time steps. The shown patterns were found by rough and preliminary parameter scans and we expect that there exist in fact many more patterns. For each example pattern, a number of questions arise: -Clusters It seems that large clusters split, leading to an intrinsic cluster size. If that is the case how is cluster size controlled and how is it determined from parameters? What determines particle density inside clusters? -Trails The observed trails appear to be a transient, unstable pattern. Is the instability due to the asymmetry in obstacle displacement and what causes it? Can other parameters produce stable trails? -Travelling bands Is there a set band wavelength and if yes, what determines it?
When are these patterns stable?
All these questions suggest that a continuous, PDE-based description of the system is crucial to understanding the observed patterns, as well as to discover others. A PDE-description has several advantages: Patterns such as travelling bands can be constructed explicitly and a stability analysis can performed. Further the PDE description is inherently an averaging process reducing the number of parameters. Lastly, since instead of numerically solving thousands of coupled ODEs, one has to solve only a few PDEs, simulations become much more efficient. The next section is therefore devoted to the derivation of the PDE-based description of the SPP-obstacle model.
large number of obstacles and SPPs

Derivation of the Macro-model
In this section, we derive a macroscopic PDE-based model for the SPPs and the obstacles. The IBM model in (1) serves as the starting point. The derivation is a twostep process: First we formally derive a kinetic description for both the SPPs and the obstacles by taking a mean-field limit. In the second step, we use a hydrodynamic scaling for the SPPs and derive equations for the SPP density and orientation. For this step, we use previous work (Degond and Motsch 2008;Degond et al. 2015). For the obstacles, we focus on a particular parameter regime and assume to have low obstacle noise and strong obstacle spring stiffness. The main technical difficulty and new derivation strategy lie in this last step. Figure 4 summarises the different derivation steps. Throughout the document, the domains of integrations are understood to mean the whole domain, unless specified otherwise.

The Mean-Field Limit
We start by defining g(x, α, t), the distribution of the SPPs at position x ∈ R n , time t ≥ 0 with direction of the self-propelled velocity α ∈ S n−1 and let f (x, y, t) be the distribution of obstacles with position x ∈ R n , tethered at y ∈ R n at time t ≥ 0.
We consider the empirical distribution associated with the dynamics of the SPPs and tethered obstacles given by system (1).
where the distributions f (x, y, t) and g(x, α, t) fulfil the following Kolmogorov-Fokker-Planck equations For the (space and time dependent) velocities, we have where we have introduced the densities of obstacles and SPPs as well as an abbreviation for densities convoluted with kernels Further f fulfils where ρ A (y) is a given, time-independent function of obstacle anchor positions.
Proof The limit is purely formal and uses standard techniques. We observe that f N and g M fulfil the equations for all N and M and then pass to the limit.

Remark 2
Note that since f and g are probabilities, they also fulfil Interpretation At this point, we have a system of coupled kinetic equations for the obstacle distribution f (x, y, t) and the SPP distribution g(x, α, t). The interactions between the obstacles and the SPPs lead to the terms of the form ∇ xρ in the speeds W and U in (6). An easy way to understand these terms is by assuming that the interaction force is of repulsive nature and purely local, in which case ∇ xρ = ∇ x ρ. We then see that the interaction force moves obstacles and SPPs in the opposite direction of the gradient of each other. The convolution with φ accounts for the potential non-locality of this interaction, which will be crucial later on. The remaining terms in W and U show the influence of the tethers and the self-propulsion for obstacles and the SPPs, respectively. In U, we also see the influence of SPP repulsion. The term involving ∇ α in (4b) reflects the effect of SPP alignment. The terms on the right-hand side of (4) are results of the stochasticity in the obstacle position (for f ) and in the SPP orientation (for g).

Scaling Assumptions
To derive the macroscopic equations for the SPP-obstacle interactions, we make a number of scaling assumptions for both the SPPs and the obstacles.
Scaling Assumptions for the SPPs Following previous work Degond and Motsch (2008), Degond et al. (2015), we introduce a small parameter ε and specify the relative order of the various terms. We mostly follow Degond et al. (2015), with a few small differences: Firstly we assume the effect of alignment to be purely local, i.e. r A = O(ε), as has been done e.g. in Degond and Motsch (2008). Alternatively one could choose a weakly non-local scaling r A = O( √ ε), which would lead to an additional viscous term in the SPP orientation equation (13b) below. As in Degond et al. (2015), we also assume the SPP self-repulsion to be purely local, i.e. r R = O(ε) and assume that However, we do not make any smallness assumption with regard to the SPP-obstacle interaction scale r I . This is because we are interested in studying the effect of the non-locality of this interaction. Otherwise we proceed as in Degond et al. (2015), i.e. assuming the alignment frequency ν and orientational diffusion d s to be of order 1/ε, and their ratio to be of order one. (6), we see that it is only the macroscopic obstacle density ρ f (x, t) that enters the SPP equation. Unfortunately, we cannot obtain a closed system for the macroscopic obstacle density ρ f (x, t) of f by integrating (4a). Instead we make assumptions about the time scales of the obstacle dynamics. From now on, we also assume to have a constant anchor density, i.e. ρ A ≡ 1 is constant in space and time. We note that the results can be generalised to non-uniform ρ A . We introduce the following quantities

Scaling Assumptions for the Obstacles From
For the derivation, we will assume both γ and δ to be small. For γ , this means that the obstacle spring relaxation time scale is small compared to the SPP domain crossing time. We will sometimes refer to this assumption as 'stiff obstacles', since it can be realised with a large spring constant κ. For δ, smallness means that the obstacle spring relaxation time scale is small compared to the obstacle diffusion time scale, which we refer to as 'low obstacle noise'. Next we rewrite (4a) as where we have defined the 'external' velocity as and the operator A y by We can rewrite the operator as where M δ (z) is a Gaussian with variance δ centred around 0, whose mass is normalised to one, i.e.
The above also shows that M δ (x − y) is in the kernel of A y .

Remark 3
Note that the rescaling of the diffusion term δ = d o γ ensures the operator A y is a Fokker-Planck-type operator. Without it, we would obtain A y ( f ) = ∇ x · [(x − y) f ], whose kernel contains Dirac deltas, making the analysis much more tedious. Eventually, however, we are interested in the small noise limit. This, of course raises several questions, which are beyond the scope of this work, e.g. does the order of the limits γ → 0 and δ → 0 matter?

The Macroscopic SPP-Obstacle Equation
Using the scaling and notation above, we now state the main result of this section, which we prove in Sect. 3.4.

where N is the von Mises-Fisher distribution defined by
Note that K d is a normalisation constant and is independent of . Further the macroscopic SPP density ρ g (x, t) and the macroscopic SPP orientation g (x, t) fulfil The constants c 1 > 0 and c 2 > 0 depend only on d = d s /ν and are defined as in Degond et al. (2015). The macroscopic obstacle density ρ f (x, t) is given by where the nonlinear term N is defined by where H(ρ g ) denotes the Hessian of the functionρ g , i.e. {H(ρ g )} i, j = ∂ i ∂ jρg , and given two n by n matrices A and B, their scalar product is defined as Equations (13a) and (13b) give the evolution for the particle density ρ g and mean orientation g , respectively. Without the term ∇ xρ f appearing in U and V in Eq. (13c), these equations correspond to the so-called Self-Organised Hydrodynamics with Repulsion (SOHR) and their derivation can be found in Degond et al. (2015). The additional terms in Eq. (13c) account for the influence of the obstacles density ρ f .
The equation for the obstacle density, expanded in the small variables δ and γ , is given in (14). It is important to note that the obstacle density given in (14) can in principle become negative, which is not physically meaningful. This is a consequence of the assumption that γ is small and indicates that the validity of the model will be limited to certain parameter regimes. We see that for infinitely strong springs, i.e. γ → 0, ρ f (x, t) ≡ ρ A ≡ 1, i.e. obstacles remain exactly at their anchor points and since those are assumed to be uniformly distributed, the obstacles have no effect on the SPPs (∇ xρ f ≡ 0). For small, but finite γ the feedback from the SPPs leads to non-uniform obstacles.

Influence of Obstacle Noise
The influence of the obstacle noise δ is contained in the order γ term in (14). We note that We see that the noise adds an additional form of non-locality. Whether the obstacle density is reduced or increased depends on whetherρ g , the convoluted SPP density at x is smaller or larger than the 'blurred', convoluted SPP densityρ g , where the amount of blurring depends on the obstacle noise. In the absence of obstacle noise, simplifies to

SPP Dynamics Deform Obstacle Volume Elements
In the absence of obstacle noise, we can rewrite (15) as where J Y is the Jacobian of the map The map Y can be interpreted as an estimate of the anchor position of an obstacle at position x moved under the influence of the SPP density. Then the determinant of the Jacobian reflects the deformation of a volume element of obstacles due to the SPPs. Note that for n = 3 det J Y contains also order γ 3 terms, for n = 2 only order γ 2 terms and lower.
Higher-Order Terms Account for SPP Movement Finally, we comment on the time derivative appearing in (14). The time derivative leads to a form of delay, i.e. the obstacles retain a memory of where SPPs were. This can be seen by Taylor expanding the SPP density in time using the time scale of obstacle relaxation γ . Then the linear terms in (15) can be written as Finally in preparation for the analytical and numerical investigation of Sects. 4 and 5, we state the following: Corollary 1 (1D equations.) Let the assumptions of Theorem. 1 hold. Then for n = 1, the equations for the SPP density ρ g (x, t) and the obstacle density ρ f (x, t) with x ∈ R and t ≥ 0 are given by where we have assumed all particles move to the right. The obstacle density up to order γ 2 is given by For δ → 0 and using only terms up to order γ , (18) simplifies to

Proof of Theorem 1
For the coarse-graining of the kinetic SPP equation (4b), we refer to previous work (Degond and Motsch 2008;Degond et al. 2015). We note that the obstacle density enters the SPP equation solely through its macroscopic density ρ f (x, t) via the interaction operator ∇ xρ f , which has a structure analogous to the SPP self-repulsion term, hence analogous techniques can be applied.
To derive an expression for the obstacle density ρ f (x, t), we formulate and prove the following Theorem: Theorem 2 Let ρ A ≡ 1 and f (x, y, t) fulfil (9) with A y ( f ) defined in (11). Let γ 1 and expand f (x, y, t) as Then the macroscopic densities defined by satisfy as δ → 0. We use the notation div = ∇ x · for the divergence of a vector field.
Proof In the following, we drop the t-dependence of most terms to increase readability. We obtain the following equations for the three highest orders of γ Let us note that (23a), (23b), and (23c) can be recast as follows: Given a function h find ψ (in a suitable functional space) such that Due to the conservation of mass property of A y , i.e. A y dx = 0, a necessary condition to warranty the existence of a solution of (24) is h dx = 0. It can be shown that the operator A y has compact resolvent on a suitable functional space and its kernel is generated by M δ (x − y), given in (12). The most important properties of the Gaussian M δ , that we will use repeatedly are Hence, we can obtain a complete characterisation of the solutions of (24) via the Fredholm alternative, namely, for any function h such that h dx = 0 there exists a unique solution ψ up to an element of the kernel of A y . For a proof of this result, consult (Aceves-Sanchez et al. 2019). Let us start by considering (23a), we search for a solution f 0 such that f 0 (x, y) dx = 1; hence, according to the results obtained for (24), the unique solution is given as where M δ is defined in (12). For the remaining two equations, we require the following scaling condition to hold, which ensures that the average mass is one, Step 1: Rescaling Next we define the functions h 1 (σ, y, t) and h 2 (σ, y, t) as This turns (23b) and (23c) into equations for h 1 (σ, y, t) and h 2 (σ, y, t). Defining B as the operator we obtain, after tedious but straightforward computations, the following relationships There are several advantages to this scaling: Firstly, the operator B is the generator of the Ornstein-Uhlenbeck stochastic process (a consequence of using σ = (x − y)/ √ δ) and we can use its well-known properties directly without having to scale by δ. Secondly, we have removed the Gaussian M δ from the equation (it cancelled). Finally, additionally scaling f 1 and f 2 by 1/ √ δ and 1/δ, respectively, turns out to be the correct choice when calculating the densities. Before we proceed to the next step, we need to collect a number of properties of B, all of which are well known and stated in "Appendix A.2".
This yields as equations for h 0 1 , h 1 1 and h 2 Note that we have used the Einstein's summation convention and that nowṽ and its derivatives are all evaluated at (y, t). Here partial derivatives are understood to act on the spatial variable, i.e. ∂ iṽ := ∂ ∂ y iṽ (y, t). The advantage of this procedure is the following: Now the right-hand sides are low-order polynomials in σ and since B only acts on σ , the equations can be solved explicitly by rewriting the right-hand sides in terms of the Hermite basis and using P2 of Lemma 4 in "Appendix A.2". This procedure yields the explicit solutions where H are the tensor Hermite polynomials defined in Lemma 4 in "Appendix A.2". Note thatṽ and all its derivatives are evaluated at (y, t) and H at σ . As equations for h 0 2 and h 1 2 , we obtain As aboveṽ and its derivatives are all evaluated at (y, t). Using the solutions for h 0 1 , h 1 1 and h 2 1 given in (29), we can solve the equations for h 0 2 and h 1 2 in the same fashion, yielding the explicit expressions Note that the solutions fulfil the scaling condition (26) since it holds that Step 3: Calculating the macroscopic moments of the obstacle density With the preparation of the two steps above, the calculation of the obstacle densities and consequently its contribution to the SPP equation becomes relatively simple. The procedure and calculations are described in "Appendix A.3". This yields (22) as claimed.
Explicit Solution for f 1 The above outlined procedure works for any given external velocityṽ(x, t), i.e. it allows to include other influences as well. For example, in future work we plan to use the derivation strategy to include the description of a fluid in which the obstacles and SPPs are immersed in. However, for this model, we can use the fact thatṽ(x, t) is in fact a conservative vector field. This allows to solve the first-order equation (23b) for f 1 (x, y, t) directly. This is covered in the following lemma, where the t dependence has been suppressed for notational convenience.

Lemma 2 Letṽ(x) be a conservative vector field, i.e. there exists a scalar function V (x)
, such that ∇ x V =ṽ, then we can write the solution to (23b) as Proof By direct calculation, we see that which shows that f 1 is indeed a solution to (23b). Finally, we have to verify the normalisation condition (26) which finishes the proof.
The above Lemma is applicable for this model of SPP-obstacle interactions since we have thatṽ i.e. we can use Lemma 2 with V (x, t) = − 1 ηρ g (x, t). We consequently find From this, we can calculate Remark 4 Note that since we see that this is consistent with (22), but contains more information about the O(δ 2 ) term.
The Macroscopic Obstacle Density Collecting the results of Theorem 2 and Lemma 2 and using the definition ofṽ given in (10), we find that the maximum order of approximation of the obstacle density we can now write is given in (14) as claimed. This finishes the proof of Theorem 1.

Analytical Insights from the 1D Macromodel.
In this section, we analyse the macroscopic model derived in Sect. 3 further to gain insights into the SPP-obstacle interactions. In particular, we use linear stability analysis to understand the onset of patterning and investigate how obstacles induce an effective SPP interaction.

Linear Stability Analysis
In this section, we investigate pattern formation for the SPP-obstacle model. We work in one space dimension, i.e. we focus on the SPP density ρ g (x, t) and the obstacle density ρ f (x, t) for x ∈ R or and t ≥ 0, whose dynamics are given by (17) and (18). Consider the steady-state solution ρ g (x, t) = ρ 0 > 0. Small perturbations of this solutions (called again ρ g ) then fulfil the linearised equation where ρ f is still given by (18).
The following propositions examine the growth or decay behaviour of perturbations of the constant steady state in dependence on their angular frequency and the resulting linear stability of the constant steady state. We consider the equation on the whole space x ∈ R and posed on an interval with periodic boundary conditions. Proposition 1 (Linear stability) Consider (33) coupled to (18) posed (a) on x ∈ R and (b) on x ∈ [0, 1] with periodic boundary conditions. (i) The system permits solutions of the form ρ g (x, t) = ρe ikx+αt , with ρ = 0, α ∈ C and k ∈ R (case a) or k ∈ 2π Z (case b) where α and k fulfil the following dispersion relation whereφ k is the Fourier transform (case a) or Fourier coefficient (case b) of the kernel φ, defined byφ where the integration domain is understood to be R (case a) or [0, 1] (case b).
(ii) The constant steady state ρ g (x, t) = ρ 0 is linearly stable iff where K = R (case a) or K = 2π Z (case b).
Proof We show the proof for case a, case b can be shown analogously. (i) Substituting the ansatz ρ g (x, t) = ρe ki x+αt into (33) is equivalent to applying the Fourier transform to the whole equation. We use the following properties of the Fourier transform and obtain an equation For the Fourier transform of ρ f , we obtain with α(k) given in (34) as claimed.
(ii) We note that the decay or growth behaviour is determined by the sign of the real part of α(k). Since the denominator will always be positive, it is sufficient to examine the numerator. This gives the result.
Corollary 2 Let ρ f be given only up to order γ and let δ → 0. Then the real part of α(k) in Proposition 1 becomes Interpretation We interpret the results of Proposition 1(ii) as indication under what conditions patterning is expected. We start by observing that in the absence of obstacle noise, δ → 0, the linear stability condition (35) simplifies to max k (kφ k ) 2 < μη γ .
Since 1 δ 1 − e −δk 2 ≤ k 2 , we observe that the obstacle noise δ > 0 has a stabilising effect. The constant on the right-hand side is critical for (in)stability. We see that SPP self-repulsion, strong obstacle springs and high obstacle friction stabilise the system. The order γ 2 approximation of the obstacle density leads to the additional terms in the denominator. It does not influence whether the constant steady-state destabilises; however, it decreases the growth or decay rate of the perturbations. The main determinant for pattern formation is the SPP-obstacle interaction kernel φ and the decay behaviour of its Fourier transform or coefficients. In case of purely local interactions,φ k is constant and we see that the we have destabilisation for all parameter values, since for large frequencies the real part of α will always become positive. This emphasises the importance of the non-locality of the SPP-obstacle interactions. Next we look at a specific case.

Example 1
We assume δ → 0 and further consider the obstacle density only up to order γ . We work on x ∈ [0, 1] with periodic boundary conditions. Further we let the microscopic SPP-obstacle interaction kernel φ be compactly supported on the interval [−r I , r I ] and yield a pushing force that decreases linearly with distance and is continuous at r I , i.e.
In this case, we can calculate the Fourier coefficients explicitly and obtain The function attains its maximum at k = π/r I and we hence find that if then the spatially constant steady state is linearly stable. The converse is in general not true, since π/r I will typically not be in 2π Z. In the case of destabilisation, we expect the pattern size P to be related to the maximum of α(k), given in Corollary 2. We observe that F(k) → 0 for k → ∞ and hence α(k) < 0 for k sufficiently large. This means high-frequency perturbations will be damped by the diffusion-like SPP self-repulsion term. Since α(0) = 0, there will typically be a well-defined maximum attained at some k = 2πl max with l max ∈ Z. We then expect that P defined by P = 1 l max will be a good indication of the expected pattern size. We numerically investigate whether this holds also far away from the constant steady state and for the IBM in Sect. 5.2.

Obstacle-Induced SPP Interaction
In this section, we show how properties of the interactions between SPPs and obstacles on the micro-level inform the properties on the macro level and find some interesting connections to equations for granular flow, porous media and aggregation equations. We focus on the simplest case, where we assume the obstacle noise to be zero and consider the obstacle equation only until order γ . Further we work in one space dimension where many calculations can be done explicitly. Then the system of interest for the SPP density ρ g (x, t) and the obstacle density ρ f (x, t) for x ∈ R and t ≥ 0 is given by (17) coupled to (19).

A Non-local Equation with Gradient Flow Structure
If we substitute ρ f given in (19) into the equation for ρ g given in (17), we obtain We now see that we have a nonlinear, non-local model with a gradient flow structure. These types of equations appear in a wide range of contexts ranging from granular flow, porous media and biological aggregation (Otto 2001;Topaz et al. 2006;Toscani 2000) and their properties are subject of intense study (Ambrosio et al. 2008;Benedetto et al. 1998;Carrillo et al. 2003). The term stemming from the SPP self-repulsion is often written as μρ = H (ρ), where H (ρ) = μ 2 ρ 2 is the SPP density of internal energy of the system. For the term stemming from the SPP-obstacle interaction, we can define the interaction kernel Note that while φ is the microscopic interaction potential between SPPs and obstacles, W can be interpreted as macroscopic obstacle-induced SPP interaction potential.

Bi-phasic Effect at the SPP Level
We now infer properties of W (the macro-interaction potential) from properties of φ (the micro-interaction potential). Note that φ > 0 or W > 0 indicate forces to the left and φ < 0 or W < 0 indicate forces to the right.

Lemma 3 (Obstacle-induced SPP interactions) Let φ(x) be an even potential. Then
φ is odd and we can define a function ϕ on [0, ∞) by using the convention that sign(0) = 0 and defining ϕ(0) := lim x→0 + ϕ(x). Let ϕ(0) = 0 and ϕ be continuous with bounded first derivative on [0, ∞). We further assume that ϕ has compact support on [0, r I ] for some r I > 0, and that ϕ and ϕ have constant but opposite sign on their support. Let W be defined as in (37). Then the following holds: (i) W is an even potential continuous on R and continuously differentiable on R\{0}. W has compact support on [−2r I , 2r I ]. (ii) W is an attractive potential for short distances, i.e. W (x) > 0 for x > 0, x small. (iii) W is a repulsive potential on (r I , 2r I ), i.e. W (x) < 0 for x ∈ (r I , 2r I ).
Proof (i): Since W is the convolution of two compactly supported, bounded functions, W is continuous. Using the definition of W and that φ is odd, we calculate Using (38), we calculate φ (x) = ϕ (|x|) + 2ϕ(0)δ(x), where δ is the Dirac delta. We therefore obtain The second term is continuous in x, since it is the convolution of two compactly supported functions, both bounded, in particular it is zero if evaluated at x = 0 due to symmetry. The first term is continuous on R\{0} hence the same is true for W . That W is compactly supported on [−2r I , 2r I ] is a consequence of the support of φ .
(ii): Using (39), we find that lim x→0 + W (x) = 2(ϕ(0)) 2 > 0, which together with the results of (i) shows that W (x) > 0 for small, but positive x. This shows that W is an attractive potential for small distances.
(iii): Let x ∈ (r I , 2r I ). Using (39), we find that By assumption, the product of ϕ and ϕ is negative, which shows that W is an repulsive potential at distances between r I and 2r I . This finishes the proof.
Example 2 (Micro-macro-potentials) We illustrate the results of the above Lemma with two examples of SPP-obstacle potentials. Using the notation introduced in (38), we consider for r ∈ [0, ∞) where H is the Heaviside function, r I > 0. ±C > 0 corresponding to attractive and repulsive SPP-obstacle interactions, respectively. The function ϕ 1 corresponds to the potential of Example 1, which is compactly supported and covered by Lemma 3, while ϕ 2 corresponds to a kernel without compact support. Figure 5a and b shows the resulting obstacle-induced SPP forces W 1 for ϕ 1 for C = −1 and C = 1, respectively. For ϕ 2 , we can see the bi-phasic behaviour directly by calculating showing that W 2 is an attractive potential for |x| < 2r I and repulsive otherwise. Note that for both examples the sign of C doesn't affect the shape of W .
Lemma 3 and Example 2 show that the SPP-obstacle interactions will have a shortranged attractive effect on SPP level, irrespective of whether the micro-interaction was attractive or repulsive. This can be understood intuitively, see Fig. 5: If the SPPs and obstacles repel each other, the obstacles that have been repelled by a group of SPPs will in turn repel other SPPs and therefore lead to further aggregation of the SPPs (Fig. 5a). On the other hand, if the SPPs attract the obstacles, the obstacles attracted by a group of SPPs will attract even more SPPs, again leading to an aggregation effect on the SPP level (Fig. 5b). Further Lemma 3 shows that if the SPP-obstacle interaction force (whether attractive or repulsive) is falling with distance, we see that in addition to the short-ranged attraction, we have a long-ranged repulsion at the SPP level as well. The second example in Example 2 suggests that this property is not limited to compactly supported functions and that Lemma 3 can be generalised to a bigger class to kernels.
These observations already give a good intuition to understand the phenomena observed in Sect. 2.2. For both the moving clusters and the travelling bands, the 1D equations (along the global alignment direction) would correspond to moving aggregates of SPPs. Both the moving clusters and the travelling bands seem to have controlled size, in particular we observed that a cluster that is too big is split into two. The above observations now give an explanation for the observed behaviour: The SPP-obstacle interaction leads to short-ranged SPP attraction and hence aggregation, however, due to the two sources of repulsion (SPP self-repulsion and obstacle-induced repulsion), clusters cannot grow too large. Next we perform 1D simulations to compare the macro-model with the IBM.

Numerical Results in 1D
In this section, we numerically solve the 1D macro-model for SPP-obstacle interactions and compare the results to both 1D IBM simulations and the analytical results of Sect. 4. Simulation details can be found in "Appendix A.1".

Comparing SOH and IBM Simulations
The Macro-SPP Obstacle Model Produces Travelling Clusters We simulate (17) coupled to (19) in 1D using periodic boundary conditions on x ∈ [0, 1] and the following parameter choices: η = 1, c 1 = 1, ζ = 8, γ = 2 × 10 −3 and μ = 5 × 10 −4 . We use a linear microscopic interaction force, i.e. the interaction kernel as defined Example 1 with r I = 0.18 and C = 0.25. As initial conditions, we use a perturbed uniform SPP density. Figure 6a shows that, indeed, moving clusters of SPPs develop, with stretches of zero density between them. The clusters seem to be relatively evenly spread. The corresponding obstacle density is minimal where the SPP density is maximal. After the clusters have been established, we inspect the space-time plot for one time unit Fig. 6b, which shows that they appear to be stably moving travelling waves of about speed one.
The Macro-SPP Obstacle Model Agrees with the IBM Next we compare to 1D IBM simulations of (1). Note that in 1D, we can disregard the orientation equation and assume all particles self-propel to the right. We use the same parameters as for the macro-model with N = M = 100 and a self-repulsion kernel yielding a linear force, dropping with distance of width r R = 0.02. As initial conditions, we use equally spaced anchor points and randomly positioned SPPs. Figure 6c shows the obstacles, their tether points and the SPPs at time t = 30. We calculate the corresponding SPP and obstacle densities from the particle positions. To that end we create a smoothed version of the empirical distribution defined analogous to (3), where the Dirac delta distributions have been replaced by 1D-Gaussians with variance 1 × 10 −4 . Note that choice of the variance is delicate, since it has to be small enough to be able to resolve the patterns and big enough to lead to meaningful averaging. The result is shown in Fig. 6d. A comparison between the simulated SPP and obstacle densities for the macro-model and IBM shows remarkable good agreement both qualitatively and quantitatively.
Higher-Order Approximations Lead to a Delay Effect The macro-model was simulated using an order γ approximation for the obstacles. To assess the effect of the order γ 2 terms without solving the full system, we proceed as follows: We substitute the measured IBM SPP density depicted in Fig. 6d into (18)  c Particle x-positions of the SPPs (red arrows with black dots) and obstacles (coloured circles, colour indicated displacement), y-positions are arbitrary. d Approximated and calculated IBM particle densities of the simulation in c, see text for details obstacle density as predicted by the model. We calculate both the order γ and order γ 2 approximations. For the latter, we need the time derivative of the SPP density, which we approximate by calculating the SPP density at the previous time step and using a forward finite difference approximation. The resulting densities are shown in Fig. 6d. We observe that measured and calculated obstacle densities agree remarkably well. Inspecting the inset in Fig. 6d, we see that the order γ approximation predicts the obstacle density minima to be precisely at the SPP density maxima; however, both the order γ 2 approximation and the actual measured IBM obstacle density have  (17), (19) showing the SPP (red, upper rows) and obstacle (blue, lower rows) densities, as well as their (constant) means (dashed black). Insets show initial condition, schematic depicts nature of microscopic interaction, red arrows indicate movement direction of the densities. Schematics in lower row: see Fig. 5 their local minima shifted backwards with respect to the SPP direction, yielding a better fit between the measured and calculated order γ 2 densities that those of the order γ approximation. This demonstrates that the derived obstacle equation allows to calculate the obstacle density for a given SPP density. It also shows that the higherorder approximation in γ is necessary if one wants to account for effects of SPP movement.

Testing Analytical Insights
Attractive and Repulsive Interactions Lead to the Same SPP Behaviour In the next numerical experiment, shown in Fig. 7 we use as initial condition a centrally placed Gaussian and inspect the moving steady state density for a repulsive (A) and an attractive (B) microscopic interaction force. We see that in both cases the resulting SPP density is the same, forming a travelling wave with a stable shape. This shape consists of a large cluster and two smaller clusters to its left and right. However, the obstacle density differs in the two cases: For an attractive potential we have obstacles clusters coinciding with the SPP clusters, whilst for the repulsive potential SPP clusters create regions of low obstacle density. The lower row compares this with the intuitive explanation of the previous section (see Fig. 5). Linear Stability Analysis Predicts Macro and IBM Patterns In Sect. 4.1, we performed a linear stability analysis for the 1D macro-equation. In Example 1, we determined the criteria for pattern formation and how to predict pattern size for a specific interaction potential shape. Now we compare these predictions to simulations of both the 1D macro-equations (17), (19) and the 1D IBM simulations by varying the size of the support of the interaction kernel r I . We use the same kernels and number of particles as above and the following parameters: η = 1, c 1 = 1, ζ = 8, γ = 2 × 10 −3 and μ = 6.7×10 −3 , C = 0.17. We start with a randomly perturbed constant initial density for the macro-model and regularly spaced anchors and randomly placed SPPs for the IBM. We compare the predicted number of peaks as calculated in Example 1 (and defined as the reciprocal of the pattern size) to the observed number of peaks at time t = 30. The result is shown in Fig. 8. We find that the analytical predictions of Sect. 4.1 agree very well with the macro model. The agreement with the IBM simulations is good as long as the macro-model gives physically meaningful (i.e. positive) obstacle densities (examples 1, 2, 3 in Fig. 8), but breaks down otherwise (example 4 in Fig. 8). This shows both that the macro-model can be used to gain insights into the IBM, but also that it is limited to certain parameter regimes.

Discussion
Summary In this work, we formulated an IBM model of the interaction of selfpropelled, collectively moving SPPs with elastically tethered obstacles. Despite the seemingly simplicity of the interactions, we found that the system can self-organise into a big variety of patterns, including travelling bands, (transiently stable) trails and size-controlled clusters. To investigate these patterns further, we derived macroscopic equations for the obstacle and SPP densities and the SPP orientation. The asymptotic regime of interest assumed γ to be small, i.e. fast obstacle spring relaxation (strong obstacle springs). The resulting continuum equations are nonlinear and contain a nonlocal interaction term. Linear stability analysis allowed to estimate pattern size from model parameters and showed which effects promote pattern formation (e.g. obstacle-SPP interaction) and which stabilise the unpatterned state (e.g. SPP self-repulsion). We found that, surprisingly, SPP dynamics are independent of whether obstacles and SPPs repel or attract each other. In particular, in Sect. 4 we discovered that the macroscopic SPP equation has gradient flow structure with a bi-phasic (short-range attractive, long-range repulsive) non-local obstacle-induced interaction kernel.
Biological Implications SPPs represent moving individuals such as animals, pathogens, bacteria, sperm, cancer or other cells. The elastically tethered obstacles mimic a complex environment, which acts on and reacts to SPPs, e.g. by repulsion. They represent, e.g. a fibrous network which is relatively fixed in space, but whose components can be pushed upon SPP interaction, after which they relax back to their original positions. Simulations revealed several patterns, such as moving clusters, where tightly packed groups of cells move together. This is a commonly observed phenomena, e.g. for groups of cancer cells invading a tissue. The key observation is that the model does not include any explicit SPP attraction, making the formation of high density patterns surprising. The subsequent analysis of the macromodel (derived in Theorem 1) gives a quantitative explanation for the observed behaviour: We found that both attractive and repulsive microscopic interactions between SPPs and obstacles cause a short-range attractive macroscopic effect on the SPP level, which leads to clustering. Clustering of organisms is ubiquitous in nature and is often attributed to direct attraction between the individuals. However, our results suggest that the apparent attraction could be indirect and is in fact mediated by the environment. In other words, it is possible the individuals feel no attraction towards each other, but will still form tight clusters. This finding is highly relevant to understanding cell clustering or swarm formation. The analysis in Proposition 1 allows to assess the formation and size of patterns (such as moving clusters or travelling bands), without the need to simulate. A key finding is that the interaction strength between SPPs and obstacles has to be large compared to the SPP self-repulsion and the tethering strength in order to create patterns. This means that, e.g. cell clustering, could be promoted by a more elastic environment. An important biological implication is that one can influence cell aggregation by modifying only the environment.
Future Work The 1D model examined in Sects. 4 and 5 cannot describe trail patterns and does not allow to distinguish between clusters and travelling bands. In the future, we plan to extend both the simulations and the analysis of the continuum model to two-and three-space dimensions. In particular, we expect new insights about pattern directionality from a linear stability analysis in higher dimensions. Our derivation relied heavily on the assumptions of smallness of γ . Mathematically this limitation manifests in the fact that the obstacle density can become negative, at which point the model becomes invalid. In the future, we would like to derive macroscopic models that are and remain well-posed for any parameter combination. This will require a different closure method of the kinetic equations. Our current model seems to be able to capture several of the observed phenomena at the IBM model, such as the travelling bands or the clusters, however, for example the trail formation pattern will most likely require an extension of the current techniques. We also plan further analysis of the continuum model, capitalising on its gradient-flow structure. Strong analytical results, such as energy dissipation estimates, exist for these type of equations, which suggests that it is possible, at least for certain cases, to construct steady states and assess their stability in a rigorous manner. Obvious extensions include 1D simulations of the SPPobstacle model using an O(γ 2 ) approximation of the obstacle density or including the positional noise. The current model describes interactions between SPPs and obstacles. In many instances, however, all components are immersed in a fluid. Past work has already studied how to derive and analyse SPP-fluid interactions ). There exist models for how fluid properties are affected if it contains immersed objects. A famous example is the Oldroyd-B model, describing the visco-elasticity of fluids filled with spring dumbbells (Oldroyd 1950). We plan to use our derivation strategy to derive equations for fluids filled with tethered obstacles and study how fluid properties such as viscosity are affected. An additional level of complexity we plan to tackle, is to combine all three components, the fluid, the obstacles and the SPPs. In this case a natural question appears: How big are the obstacles compared to the SPPs. The flexible techniques developed in this work will allow to answer this question by performing the coarse-graining at different levels.
Applying the findings to biological systems, such as sperm movement, will require careful parametrisation of the model. An advantage of the IBM formulation is that model parameters are relatively easy to obtain from experiments: Diffusion constant can be estimated using the Stokes-Einstein formula and interaction radii can be linked to object sizes. Other key properties, such as how forces behave with respect to distance, could be measured by observing how an individual SPP reacts to an individual obstacle. The model then allows to predict the result of the dynamics of large groups of SPPs and obstacles. Note that H j (s) are the (probabilistic) Hermite polynomials. Let us consider the following L 2 -weighted space

IBM Simulation Videos
where M 1 is defined in (12) taking δ = 1. For any two functions f and g in X , we define their weighted inner product by f , g X := R 3 f (σ )g(σ )M 1 (σ ) dσ.
We then have the following properties: We have used the notation i! = i 1 !i 2 !i 3 ! and |i| = i 1 + i 2 + i 3 . Note that P1 shows that H i are orthogonal with respect to the inner product ·, · X and P2 states that H i are eigenfunctions with eigenvalue −|i|. In the product rules P4 and P5, e i denotes the i-th unit vector in R 3 .

A.3 Calculation of the Obstacle Density
In this section, we detail the calculations of 0-th order moment of f , ρ f , in terms of expansions with respect to γ and δ. As outlined in the main text, we will perform the following steps: 1. Perform the change of variables √ δσ = x − y. This changes the integrand to be proportional to M 1 (σ )h i (σ, x − √ δσ, t) for ρ f i . 2. Next we Taylor expand h i (σ, x − √ δσ, t) around δ = 0 using the expansions of above. 3. Then we calculate the contributions using the scaling condition (31), the orthogonality of the Hermite polynomials (P1) and the product rule (P4) of Lemma 4.
The following result will be helpful for the subsequent calculations.
Lemma 5 Let h k 1 and h k 2 , for k = 0, 1, . . ., be the solutions of the above expansion and let their representations w.r.t the basis of Hermite polynomials be given by where i is a multiindex and a k i and b k i are functions of y and t. Then it holds that a k i ≡ 0 for |i| mod 2 = k mod 2, or |i| > k + 1 b k i ≡ 0 for |i| mod 2 = k mod 2, or |i| > k + 2