A multi-layer model for self-propelled disks interacting through alignment and volume exclusion

We present an individual-based model describing disk-like self-propelled particles moving inside parallel planes. The disk directions of motion follow alignment rules inside each layer. Additionally, the disks are subject to interactions with those of the neighboring layers arising from volume exclusion constraints. These interactions affect the disk inclinations with respect to the plane of motion. We formally de-rive a macroscopic model composed of planar Self-Organized Hydrodynamic (SOH) models describing the transport of mass and evolution of mean direction of motion of the disks in each plane, supplemented with transport equations for the mean disk inclination. These planar models are coupled due to the interactions with the neighboring planes. Numerical comparisons between the individual-based and macroscopic models are carried out. These models could be applicable, for instance, to describe sperm-cell collective dynamics.


Introduction
Collective motion in systems of self-propelled particles is the subject of a vast literature. How collective motion emerges from the underlying local interactions between the agents is still poorly understood. The interactions are either of cognitive nature (such as in birds, mammals) 29 and/or are mediated by a surrounding fluid (such as in swimming bacteria, sperm cells, etc.). 18 There are several competing strategies to model collective dynamics. One strategy relies on individual-based models 1,[5][6][7][8]17,19,21,26 that describe how the position and velocity of each individual evolves in the course of time. Another strategy relies on continuum models 2,3,25,27 describing the system by locally averaged quantities such as the mean density, mean velocity, etc. An intermediate category of models consist of kinetic models 15 which describe the individual motions in a probabilistic way. These different types of models can be connected one to another as kinetic models can be seen as resulting from an infinite particle number limit of individual-based models, while macroscopic models follow from a hydrodynamic or diffusive rescaling of the kinetic models and subsequently passing to the limit of a large rescaling factor. We refer to Ref. 12 for an illustration of this methodology in the case of the Vicsek model (see also below).
In this paper, we are concerned with finding a suitable modeling framework for the three-dimensional motion of spermatozoa in the seminal plasma. As there are evidence that sperm-cell motion in the most common experiments is mostly planar, 24 we propose a multi-layer model where the motion of sperm-cells is planar and sperm-cells may interact with other sperm-cells of the same layer or of the two neighboring layers. Sperm-cell concentration in raw sperm is incredibly high (as large as 5 × 10 9 cm −3 ) and the use of individual-based models to reproduce real sperm-cell experiments is intractable. It is therefore necessary to derive a macroscopic model describing the collective motion of the cells within each layer and their interactions with the neighboring layers.
Spermatozoa can be assimilated to two-dimensional discs. Indeed, as regard to the occupied volume, flagella can be neglected and spermatozoa reduced to their head. Although the heads resemble flat ellipsoids, we simply model them as infinitely thin flat discs. The acting flagellum produces almost constant propulsion. Thus, it is a good approximation to suppose that all sperm-cells move with the same constant speed and that only the velocity direction is subject to changes. Finally, each disk possesses some inclination with respect to its plane of motion, measured by an inclination angle. Therefore, the position, velocity and attitude of each disk can be described by the position of its center of mass, its velocity direction and its inclination angle.
Interactions between sperm-cells are mostly hydrodynamic interactions (due to the perturbation of the fluid velocity induced by the motion of the cells) and volume exclusion (or steric) interactions (due to the impossibility that two sperm-cells overlap). Modeling hydrodynamic and steric interactions within dense suspensions of active particles is a difficult subject. However, it has been shown 23 that for selfpropelled elongated particles, such interactions simply result in local alignment of the particles with their neighbors. Therefore, we assume that all these interactions can be lumped into a local alignment interaction with the close neighbors.
As already mentioned, we assume that particle motion is planar and takes place in parallel two-dimensional layers. Each particle belongs to one layer for all times without the possibility to change its layer. For the reasons outlined above, spermatozoa interact inside these 2D layers by alignment of both their velocity and inclination with those of their close neighbors. Specifically, we consider the time-continuous version of the Vicsek microscopic model as proposed in Ref. 12, where each particle tends to align with the mean direction of its neighbors up to a small Brownian perturbation. The original model was proposed by Vicsek et al. 28 and several variants have been proposed. 9,11,16,22 The inclination alignment dynamics follows a similar rule with the exception that the interaction is nematic (i.e. two inclination angles differing by a multiple of π lead to the same disk attitude). The combined alignment dynamics in the velocity-inclination variables thus differs from the 3D Vicsek dynamics (specifically, in the 3D Vicsek, velocity belongs to a 2Dsphere whereas here the pair (velocity, inclination) belongs to a 2D torus).
The different layers also interact via the volume-exclusion constraint. Indeed, due to the inclination of the disks, the spermatozoa of one layer exert a friction on the spermatozoa of the nearby layers. This interaction results in increasing or decreasing the inclination of the spermatozoa in these layers. A given layer thus acts on the neighboring ones in a similar way as the wind does on plant canopies. 14 The involved mechanical forces depend on the geometric configuration of the discs: it thus depends on their respective inclinations and also on their respective velocities. It results in alignment of velocities and repulsion of inclination angles. Layers are thus coupled and this coupling depends on the so-called overlap function that quantifies the distance between layers.
We then consider a mean-field kinetic version of the model. This equation provides the time evolution of the distribution function in phase space (position, velocity and inclination) and takes the form of a Fokker-Planck equation. Here, it is formally derived in the limit of an infinite number of interacting particles. Although the mathematical validity of the mean-field limit has been established for the Vicsek model, 4 it is still open for the present model and our result so far is only formal. We perform a spatio-temporal hydrodynamic rescaling of the mean-field kinetic model, considering that the intra-layer interaction scales are much smaller than those of the inter-layer interactions and that the latter occur at the same scales as the macroscopic evolution of the system. There results a singularly perturbed Fokker-Planck equation, involving a small parameter ε measuring the ratio of the small (microscopic) scale to the large (macroscopic) one.
The macroscopic description of the system is found by letting ε to zero in the singularly perturbed mean-field kinetic model. We first need to find the equilibria associated to the Fokker-Planck operator. We show that these are given by products of von Mises distributions in the velocity and inclination angles respectively. von Mises distributions are the natural analog of Gaussian distributions for probabilities on the sphere. We then need to integrate the equation against the collisional invariants. However, as noticed in Ref. 12, only mass is a collisional invariant and we are thus lacking two more collisional invariants to obtain the dynamics on the velocity and inclination. Following Ref. 12, we have to introduce the "generalized collisional invariants" of this operator: these are collisional invariants valid only on functions with prescribed mean velocity and mean inclination. We are then able to derive the macroscopic model.
The obtained model consists of a continuity equation for the density ρ and two evolution equations for the mean velocity angleφ ∈ [0, 2π] and mean inclination angleθ ∈ [0, π]. All these quantities are indexed by h corresponding to the hth layer. The model is written as: where V (φ) = (cosφ, sinφ) T denotes the velocity vector. The constants will be defined further. The left-hand sides of (1.1)-(1.2) form the SOH (self-organized hydrodynamics) model describing the Vicsek dynamics at the macroscopic level. 12 They respectively account for the conservation of mass and convection of the mean velocity angle. The left-hand side of (1.3) describes the advection of the inclination with the same advection velocity as for the mass. Finally, the right-hand sides of (1.2)-(1.3) describe the inter-layer interactions. The obtained model is thus a hyperbolic system (like the SOH model) with source terms that couple the velocity and inclination dynamics. We remark that, once the velocities of all the different layers are colinear, the source terms vanish. The model thus simplifies into a superposition of standard SOH models in each layer. However, before reaching an equilibrium, the interplay between inclination and velocity crucially determines which equilibrium velocities and inclinations will be attained.
To validate the macroscopic model, numerical simulations are performed and compared with those of the particle model. With this aim, we adapt the numerical relaxation method 20 designed for the Vicsek model. The numerical simulations show that the macroscopic model captures the velocity alignment between the different layers quite well. However, some differences in the inclination dynamics appear. Indeed, some transient "meta-stable" configurations arise during the course of time. They are more rapidly left away by the microscopic dynamics than by the macroscopic ones, probably because of the stochastic fluctuations associated with the microscopic dynamics. This effect due to the finiteness of the particle number in the microscopic dynamics could probably be reproduced by including a stochastic term in the macroscopic model such as the one proposed in Ref. 27.
The outline of this paper is as follows. In Sec. 2, we introduce the microscopic model. In Sec. 3, we present the mean-field limit and the hydrodynamic scaling. In Sec. 4, we provide the derivation of the macroscopic model. In Sec. 5, we compare the microscopic and macroscopic dynamics on several test-cases. A discussion of the results is provided in Sec. 6. Appendix A provides complements on nematic alignment modeling, Appendix B gives some properties of the coefficients of the macroscopic model and Appendix C describes the numerical schemes used for the simulations.

Microscopic Model
Spermatozoa are represented by discs of radius R moving in different layers. In the first approximation, the layers, indexed by Z, are copies of the R 2 plane: layer h denotes the plane {(x, y, z) ∈ R 3 , z = hd} where d > 0 denotes the inter-distance between the layers. With no interactions between layers, spermatozoa are supposed to be orthogonal to the layer and follow the Vicsek dynamics. However, if the layer inter-distance is lower than the disc radius R, spermatozoa of layer h will exert a force on the spermatozoa of the neighboring layers h − 1 and h + 1: they may force them to incline (with respect to the layer plane).
We consider N discs in R 3 labeled by k ∈ {1, . . . , N}: each disc is contained into one layer and thus disc k has a permanent altitude h k ∈ Z. The two-dimensional movement into the layer is described by the position of its center of mass X k (t) ∈ R 2 and the velocity orientation of its center of mass V k (t) ∈ S 1 : we indeed consider that all particles move at the same speed c > 0. The velocity of the particle is thus given by cV k . We introduce the angle ϕ k (t) of V k (t) with respect to a reference axis, so that V k (t) = V (ϕ k (t)) with V (ϕ) = (cos ϕ, sin ϕ).
Concerning the configuration of the disc in space, we suppose that the disc moves in one direction contained in its plane: the orientation V k (t) ∈ S 1 belongs to the disc plane. The angle of this plane with respect to the z-axis is denoted θ k (t) ∈ R/[0, π]. An angle θ k = 0 or π means that the disc is perpendicular to the plane while an angle θ k = ±π/2 means that the whole disc lies in the layer plane.

Dynamics for the centers of masses
The centers of masses follow a Vicsek-like dynamics as introduced in Ref. 12. The dynamics of the positions and the velocities are given by the following equations: Two dynamics are in competition: alignment and diffusion. Each particle of a given layer h k tends to align with a direction V (φ tot k (t)) with an intensity ν supposed constant. The direction V (φ tot k (t)) is defined as a weighted mean direction of the neighbors particle k in layers h k , h k ± 1 within the disc of radius R 1 : where J ϕ k (t) denotes the contribution of neighbors belonging to the same layer h k , J ϕ,w,nb k (t) denotes the contribution of neighbors belonging to layers h k ± 1 and β quantifies their relative involvement. Superscripts "nb" means neighboring layers and "w" means weighted. Indeed, due to steric constraints within layers, directions of neighbors of layers h k ± 1 are weighted according to their inclination. Supposing distance h between layers is smaller than twice the particle radius, h < 2R, we define the overlap function g by: This function quantifies the overlapping area of two discs in theẑ-direction. The particle direction is also submitted to a Brownian motion B ϕ,k t with diffusion coefficient D > 0. We here neglect congestion forces and we also neglect alignment between discs of different layers.

Dynamics of the disk orientations
The disc angle dynamics follows the following torque balance equation in the overdamped regime a : a The torque balance equation reads: I d 2 θ k dt 2 +C dθ k dt +K sin(2(θ k −θ k )) = T k , where I is the moment of inertia of the disc with respect to its longitudinal axis (parallel to V k ) and C is the dissipation coefficient. In the over-damped regime, I is negligible and we recover (2.6) with K =K/C. where the alignment dynamics inside each layer is in competition with steric forces between layers. The factor 2 inside the sine function takes into account that the inclination interaction between disks is a nematic one, i.e. orientations θ k =θ k or θ k = −θ k are equivalent. In this way, Eq. (2.6) preserves the fact that θ k is defined modulo π. K is the rotational stiffness and δ > 0 is a diffusion coefficient. The inclination of each particle in layer h k tends to align with the mean inclinationθ k of neighboring particles of the same layer in the disc of radius R 2 . The mean inclination is defined through the following average, which corresponds to a nematic alignment mechanism (see Appendix A): Particles on neighboring layers h k ± 1 exert a steric force on particles of layer h k . T k (t) is the sum of the weighted torques T kj (t) (with respect to the longitudinal axis) of the forces exerted by the discs: whereẑ is the unit vector in the vertical direction and g(θ j , θ k ) is the weight defined in (2.5). Indeed, the force F kj exerted by disc j on disc k is supposed to be the projection of the velocity direction V j on the orthogonal plane to V k : where µ is a mobility coefficient. Then the torque T kj of the force F kj with respect to the longitudinal axis of the disc is given by: and using the approximation

Mean-field kinetic model
We introduce the distribution function in phase space: f (x, ϕ, θ, h, t), where f is 2πperiodic in ϕ and π-periodic in θ. Note that h is still a discrete parameter numbering We can easily show the following expansion: where the last quantities are the localized mean inclination angle and momentum. Therefore, system (3.1)-(3.5) becomes, dropping O(ε 2 ) terms: where We suppose also that the interaction between the layers happens on large time scale. Therefore, we write: 1 ε µRπR 2 3 = µ = O (1). Inserting this ansatz, we obtain: Moreover, we suppose that β = εβ with β = O(1). We thus have the following expansion b : and Taking the cross product of the previous expansion with V (ϕ), we easily obtain: Therefore, regrouping the O(ε) terms together, Eq. (3.9) becomes: System (3.10)-(3.16) is the starting point for the derivation of the macroscopic model.
b For any vectors a, b ∈ R 2 and ε > 0, we have:

Equilibria
We want to take the limit ε → 0 in system (3.10)- (3.16). Therefore, assuming that the distribution function f ε converge to a limit denoted by f 0 as ε → 0, this limit satisfies the equilibrium condition Q(f 0 ) = 0 where We define the von Mises-Fisher (VMF) distribution with periodicity 2π/n, n ∈ N\{0}, concentration parameter κ > 0 and direction ψ 0 ∈ R/(2πnZ) by: We introduce We show the following: We have
We define We have: (4.7) In particular, we havē This proposition shows that for the distribution ρMφ ,θ , ρ is its local density andφ,θ are its mean direction of motion and mean inclination respectively.
is given by: Now, let f be an equilibrium, i.e. Q(f ) = 0. Then, This shows that there exists Reciprocally, we show that f = ρMφ ,θ are such that Q(f ) = 0. Indeed, this follows from (4.8) and (4.9) if we show that for f = ρMφ ,θ ,θ f =θ andφ f =φ. But this follows in turn from Proposition 4.2.
According to this proposition, there exist Now, the goal is to find equations for (ρ,φ,θ) as functions of (x, h, t).

Collisional invariants
Equation (3.15) can be written as In the limit ε → 0, the right-hand side of (4.11) is singular. The goal of a collision invariant is to cancel this singular term through integration in (ϕ, θ) against a suitable test function. For this purpose, we define: We denote by C the space of CI. It is a vector space.
Here, clearly, C contains the constant functions, since, as seen from (4.3), we have Q(f )dϕ dθ = 0. However, no other CI appears obviously from (4.3). The use of the constant functions as CI already gives the mass conservation equation. Indeed, integrating (4.11) with respect to (ϕ, θ), we get In this equation, the 1/ε singularity has disappeared and the limit ε → 0 of Eq. (4.12) leads to an equation for f 0 . However, as seen from (4.10), f 0 depends on three scalar quantities: ρ,θ andφ but (4.12) is only one single scalar equation. Therefore, it is not sufficient to determine the dynamics of f 0 . For this reason, we look for a weaker invariant concept, that of generalized collision invariant, as defined in the next section.

Generalized collisional invariants
For a given pair (φ,θ) ∈ R/2πZ × R/πZ, we define the operator Q(φ,θ; f ) by . (4.13) We note that f → Q(φ,θ; f ) is a linear operator and that (4.14) We define the generalized collisional invariants by the following: We denote by G the space of GCI associated with (φ,θ). It is a vector space.

Macroscopic equations
We obtain the macroscopic dynamics by integrating system (3.9) against the GCI. The resulting equations are presented in the following proposition:

Macroscopic equilibria
One simple macroscopic equilibrium consists in layers with the same vector velocity fields (or opposite vector field). In that case, inclinations have no impact on the dynamics and are simply transported by the velocity flow. In particular, we have: define a homogeneous macroscopic equilibria.
The case of opposite vector flows may be unstable since a small deviation from the equilibria leads to the global alignment of the layers. In particular, numerical simulations (see Sec. 5.1.2) capture only equilibria with the same velocity in each layer. The question whether other stable macroscopic equilibria exist remains open.

Numerical Experiments
In this section, we compare numerical simulations of both the microscopic and macroscopic models. The numerical methods are variations of those presented in Ref. 20: the microscopic model is solved with an implicit scheme and the macroscopic model is solved using the splitting method. The two methods are detailed in Appendix C.

Convergence to equilibria
We consider three layers and each layer contains 600 particles. The particle positions are uniformly randomly distributed on the square [0, L x ] × [0, L y ] with L x = L y = 1, the particle velocity direction angles ϕ are uniformly randomly distributed on R/[0, 2π] and the particle inclinations θ are uniformly randomly distributed on R/[0, π]. We first consider a homogeneous test-case: the interaction radii R 1 , R 2 and R 3 all equal L x /2 and each particle thus interacts with (almost) all the others.
We first consider non-interacting layers supposing h > 2R. In Fig. 2, we plot the distributions of ϕ and θ of the three layers. We also plot the von Mises distributions: whereφ andθ are the mean velocity angle and mean inclination angle of the particle of the th layer: Both velocity and inclination distribution are in good agreement with the von Mises distributions. In Fig. 3, we present the time evolution of the mean anglesφ and θ . Since there are no layer interactions, these mean angles are almost constant in time up to stochastic fluctuations.

Interactions between layers
We still consider three layers but the number of particles per layer equals 25,000. The particle positions are uniformly randomly distributed on the square [0, L x ] × [0, L y ] with L x = L y = 1 and the interaction radii R 1 , R 2 , R 3 equal 0.02: there are  in average 31 neighboring particles. We include layer interactions: we suppose that the inter-layer distance h equals the particle radius R = 0.02. The layer interaction coefficients are chosen as follows: µ = 3 and ν = 0.5. Particle velocity and inclination angles are randomly distributed according to their respective von Mises distribution. Initial mean velocity angle and mean inclination angle for the three layers are chosen as follows:φ We also perform a time rescaling in the microscopic model. Let ε > 0 and consider the following microscopic parameters: with ε = 0.1. Consequently, we choose a macroscopic time scale. We compare particle simulations with macroscopic simulation. For the macroscopic model, we thus consider a constant initial density in each layer given by: and the uniform initial values ofφ andθ given by (5.1) and (5.2). The macroscopic layer interaction coefficients are given by: with no ε, since the particle simulation already considers the macroscopic time scale. Figure 4(top) depicts the time evolution of the mean velocity and mean inclination for both the microscopic simulation (dashed line) and macroscopic simulation (continuous line). Microscopic simulations are averaged of 20-particle simulations. Let us first describe the macroscopic dynamics. The dynamics can be split into two steps: during the first step, up to time t ≈ 2, layer 2 mainly interacts with layer 1 since the overlap function is more important between these two layers. This leads to a first relaxation dynamics that make the mean velocities of these two layers align. Then, during the second step, after time t = 2, interactions between layers 2 and 3 become predominant and a second relaxation dynamics occur that leads to alignment of the three layers. We then note that microscopic and macroscopic simulations coincide during the first 1.5 time unit (including the first relaxation mechanism). This confirms that the macroscopic model captures the right interaction time scale. However, we see that, due to the finite number of particles, stochastic fluctuations make the second relaxation occur earlier around time t = 3 (micro) instead of t = 3.5 (macro). Looking at Fig. 4(bottom), where 10-particle simulations are plotted and compared to macroscopic simulation, we see that the time of the second relaxation depends on the simulation and always occurs before the macroscopic relaxation. As noted in Sec. 4.5, once the particle velocities of the three layers are aligned, homogeneous inclination angles per layer define equilibria. Therefore, the time of the second relaxation step strongly determines the final inclinations. This explains the large deviation between macro and micro simulations after time t ≥ 4 as regard to inclinations.
We now conserve the same parameters except that h = 0.0205 > R = 0.02 and, consequently, some inclination configuration prevent layers from interacting. In Fig. 5, we plot the time evolution of the mean velocity angle and mean inclination angle. We observe that, in the macroscopic simulation, a slight increase of the interlayer distance results in large time translation of the second relaxation step going from t = 3.5 (Fig. 4) to t = 7.5 (Fig. 5). This highlights the meta-stability of the system between the two relaxation steps. Concerning the particle simulations (in dashed lines), the slight increase of the inter-layer distance does not result in a so much increase of the second relaxation time (it goes from t = 3 (Fig. 4) to only t = 3.5 (Fig. 5)). The second relaxation thus occurs two times earlier than predicted by the macroscopic simulation. Indeed, due to stochastic fluctuations, some particles interact instead of remaining in non-interacting configuration. Therefore, this is the stochastic fluctuations that impact the long-term dynamics of the model.

Inhomogeneous simulations
We now consider three layers on the square domain [0, L x ] × [0, L y ] with L x = L y = 10. As in Ref. 9, we are interested in Taylor-Green vortex initial condition. Initial densities are taken uniform equal to 10 3 . Velocity angles are given by: where w denotes the angle between vectors w ∈ R 2 and (1, 0) T , (u, v) is the vector defined by: u(x, y) = 1 3 sin π 5 x cos π 5 y + 1 3 sin 3π 10 x cos 3π 10 y and (x ,ȳ ) are translation of (x, y): with translation vectors (t x , t y ) given by: As regard to the velocity initial condition, each layer is thus the translation of a normalized Taylor-Green vortex. Consequently, layer velocity fields are not initially the same and alignment dynamics should arise. Finally, inclination angles are taken uniform with the same value as the previous test-case (5.2).

Non-interacting layers
We first consider non-interacting layers: h > 2R. This test-case thus reduces to a simulation of the SOH model. In Fig. 6, we plot the space distribution of density, velocity and inclination for both the particle (left figures) and macroscopic simulations (right figures) for layer 2 (first and third layers are identical up to translation). For the particle simulation, we consider that each layer contains 10 5 particles. The interaction radii R 1 , R 2 and R 3 equal 0.04. Therefore, the particles have in averaged πR 2 1 10 3 ≈ 5 neighboring particles in each layer at the beginning of the simulation. We consider the same rescaling (5.3) with ε = 0.1. Densities are computed on a grid with space steps equal to ∆x = ∆y = 0.2 and as regard to the particle simulation, they are averaged over 20 runs of the test-case. We use the same time step for both macroscopic and microscopic simulations.
In Fig. 6(top), we observe clustering for both microscopic and macroscopic simulations in region of negative divergence flow. The density in layer 2 has maximal value equal to ρ 2 ∞ = 6520 for the microscopic simulations and ρ 2 ∞ = 7547 for the macroscopic simulation. Consequently, the number of neighboring particles  is multiplied by a factor 10 in these regions. We also observe that the vortex are better conserved with the particle simulations. Concerning the velocity field, the microscopic velocity is obtained by dividing the local momentum by the local density. The velocity vectors should be of size c 1 but this is not the case in low density regions due to the small number of particles. This partly explains the observed differences between the microscopic and macroscopic simulations. However, there is quite a good overall agreement between the two simulations. Figure 6(bottom) shows the cosine (in absolute value) of the inclination angle: as layer interactions do not occur, it remains uniform.

Interacting layers
We then consider interacting layers when setting h = R. The layer interaction parameters are taken equal to: β = 2, µ = 20. The other parameters are the same as previously. The parameter µ is chosen in such a way that the interaction coefficient µRπR 2 3 ρ 0 ≈ 0.5 is of the same order as in Sec. 5.1.2. We compare the average of 20-particle simulations with one macroscopic simulation on Figs. 7 and 8. Figure 7 depicts the density and the velocity vector field for the three superposed layers. As in the previous test-case, we observe that the microscopic and macroscopic simulations are in good agreement. The velocity vector fields of the three layers are mostly aligned and consequently the density concentrations are localized almost in the same regions. This is particularly true for the macroscopic simulations (right). Concerning the microscopic simulations (left), we still observe some differences between the layers.
In Fig. 8, we represent the cosine of the inclination angle (in absolute value). The inclination for particle simulations is obtained by computing the local mean inclination angle with formula (2.7). Due to layer interactions, the inclination is no more uniform. Contrary to the velocity vector field, we observe large differences between the macroscopic and particle inclinations. This could be a consequence of the differences pointed out in Sec. 5.1.2. Note that regions with aligned inclinations do not necessarily match regions of uniform densities. Finally, the inter-layer interactions on inclinations affect in return the density and velocity vector field. This is particularly clear when comparing the density of layer 2 with the one of the non-interacting test-case ( Fig. 6(top)) for the macroscopic simulation. All the symmetries inherited from the initial distribution have been diluted by the inclination/velocity inter-layer interactions.

Conclusion and Discussion
In this paper, we have proposed an individual-based model of self-propelled disklike particles interacting through alignment and volume exclusion. This model is intended to provide a framework for modeling collective sperm-cell dynamics. Particle motion is supposed to be confined in two-dimensional planar layers. Particle interactions between nearby layers contribute to modify the disk inclinations, which generates a coupling between inclinations and motion. We have then derived a continuum model from this individual-based model. It describes the evolution of the local density, mean velocity direction and mean inclination of the disks in the various layers. Numerical simulations have shown a good agreement between the continuum model and the individual-based one, but has also highlighted some differences.
There are many possible directions to expand the current work and make it more realistic. At the individual-based level, the description of the agents and the interaction rules could be improved for a better account of actual sperm-cell motion. For instance, the assumption that all sperm-cells have the same constant velocity is obviously unrealistic. In any semen sample, there is always a certain proportion of dead sperm-cells and of less motile ones. This could be accounted for by allowing the particle velocities to span a certain range of values. The shape of the head could be improved from the current infinitely thin disk to finite thickness ellipsoids. The inclination interaction could involve a density dependency as it is more difficult to fit actual disks in one layer if they are inclined towards the plane than if they stand vertically. Finally, one could also imagine a process by which particles would change layers. At the level of the continuum model, one major improvement should be to add a random fluctuation term in order to account for finite system size effects, similar to Ref. 27. We believe that adding such a term would help achieve a better match between the continuum model and the individual-based one. Other improvements would consist in adding a spatial diffusion to retain some of the nonlocality of the alignment interaction, similar to Ref. 11 or adding a layer-changing term. Finally, for both the individual-based and continuum models, the model parameters should be calibrated by close comparisons with biological data.

Appendix A. Nematic Alignment
Nematic alignment consists in alignment with the mean direction. The mean direction can be defined as the eigenspace of the averaged projection matrix: Thenθ is defined modulo π/2. Therefore, the orthogonal vectors (cosθ, sinθ) T and (cos(θ + π/2), sin(θ + π/2)) T are both eigenvectors. Using the relation θ is now defined modulo π. We recover the definition of Eq. (2.7).

Appendix B. Coefficients of the Macroscopic Model
From Ref. 12, coefficients c 1 and c 2 are positive and bounded. We have the following expansion (in the limit κ 1 → 0, that corresponds to large diffusion compared to alignment): The following proposition asserts similar results for coefficients c 3 and c 4 .
Proposition B.1. Constants c 3 and c 4 can be written: π 0 e κ1 cos udu π 0 e −κ1 cos udu , π 0 e κ2 cos 2udu π 0 e −κ2 cos 2udu . In particular, they are positive and lower than 1. We have the following Taylor expansions, as κ 1 → 0 and κ 2 → 0: Proof. By integration by part: Then, using expression (4.23), we easily get the expected expression for c 3 . The same manipulations lead to the expression for c 4 . The positivity of c 3 and c 4 then results from the Cauchy-Schwarz inequality.
Moreover, the weight functions gM 2 M 2 and gM 2 M 2 ∂ θ I 2 are also bounded.
In the last expression, both terms have absolute value lower than 1.

Appendix C. Numerical Schemes
For the sake of completeness, we here recall the numerical scheme used for the numerical simulations.