Continuum model for linked fibers with alignment interactions

We introduce an individual-based model for fiber elements having the ability to cross-link or unlink each other and to align with each other at the cross links. We first formally derive a kinetic model for the fiber and cross-links distribution functions. We then consider the fast linking/unlinking regime in which the model can be reduced to the fiber distribution function only and investigate its diffusion limit. The resulting macroscopic model consists of a system of nonlinear diffusion equations for the fiber density and mean orientation. In the case of a homogeneous fiber density, we show that the model is elliptic.


Introduction
The topic of complex systems is attracting an increasingly abundant literature, due to its paramount importance in life and social sciences. Complex systems consist of a large number of agents interacting through local interactions only and yet able to self-organize into large-scale coherent structures and collective motion [36]. Among examples of interactions leading to collective motion, the alignment interaction has been the subject of many studies since the seminal work of Vicsek and co-authors [35]. In Vicsek's model, self-propelled point particles tend to align with their neighbors up to some noise. Vicsek's particles are polar: they carry a definite direction and orientation defined by the unit vector of their propulsion velocity. Their alignment interaction is also polar in the sense that a particle moving in an opposite direction to its neighbors will eventually reverse its direction of motion. However, other alignment rules have been studied as well. Polar particles can be subjected to nematic alignment. In this case, a particle moving in an opposite direction to its neighbors will not reverse its direction of motion, as opposed to the polar alignment case. Nematic alignment has been used as a model for the volume exclusion interaction [5,20,29] .Particles can also be apolar, for instance if they randomly reverse their direction of motion. Apolar particles interacting through nematic alignment have been proposed as a model for vibrating rods [6], or fiber networks [1]. In the related field of nematic liquid crystals, volume exclusion interactions between rod-like particles are also modelled as an alignment force [18,24,28]. But additionally, the molecules are convected by the background solvent and are subjected to rotation by the fluid shear. Additionally, they contribute to the fluid dynamics of the liquid solvent through an additional extra-stress tensor. Usually, the polymer chains are supposed of fixed length, although lately, models of polymer chains of variables lengths have appeared [12].
In the present work, we are interested in a system consisting of fibers (or polymer chains) of variable lengths. This model aims to describe the network of collagen fibers in a fibrous tissue. We model fiber length variation (through polymerization / depolymerization) as well as the ability for the fibers to establish cross-links between them by the same basic rules described as follows. We assume the existence of a fiber unit element (or monomer) modeled as a line segment of fixed length L. We suppose that two fiber elements that cross each-other may form a link, thereby creating a longer fiber. There is no limit to the number of cross-links a given fiber can make. Therefore, the fibers have the ability to branch off and to achieve complex network topologies. We include fiber resistance to bending by assuming the existence of torque which, in the absence of any other force, makes the two linked fiber elements align with each other. Fibers are also subject to random positional and orientational noise and to external positional and orientational potential forces. Finally, cross-links may also be removed to model possible fiber breakage or depolymerization.
Our model features apolar fiber particles (since they are not self-propelled), interacting through nematic alignment with the other fibers they are linked to. Thus, the model bears analogies with previous models of apolar particles interacting through nematic alignment [6,1]. However, the interaction network topology (which keeps track of which fiber pairs are cross-linked) is different, as ours is determined by the distribution of cross-links. The fact that this network topology changes with time through dynamic cross-linking or unlinking processes is one specific feature of the present work. In the absence of crosslink remodeling, i.e. when the cross-links lifetime is infinite and no new cross-links is created, each connected component of the fiber network can be seen as an unstretchable elastic string since all connected fiber elements will spontaneously align with each other. However, cross-link removal or creation events (supposed to occur at Poisson distributed random times) introduce a fluid-like component to the rheology of the fibers, thereby confering some visco-elastic character to the medium. Cross-link-governed statics and dynamics of fiber networks have been intensely studied in the literature [3,8,9,21,27] . However, most models consider passive cross-links which only act on the fibers by a spring-like attractive force. Here, our description introduces active links which tend to align the two fibers with each other. By doing so, we are also able to take into account fiber breakage, elongation and branching just in addition to and in the same way as fiber linking/unlinking because cross-linked fiber elements can be seen as two parts of the same fiber. Another difference from previous literature is that fibers in our model are subject to noise making the system more akin to a fluid or a gas than to a solid. By contrast to classical polymeric fluid studies, we do not assume that the fibers are transported by a fluid and modify its rheological properties but this feature could be added in future work.
This model was first introduced in Ref. [30] where it was coupled with the dynamics of spherical particles modelling cells. This model has been built to describe the selforganization of the adipose tissue, where spheres represent adipocytes and fibers, the surrounding collagen fibers. In this work, we demonstrated that the interaction between cells and fibers led to the spontaneous formation of cell clusters of ovoid shape akin to the adipose lobules that form the functional subunits of the adipose tissue. In Ref. [30], only a discrete Individual-Based Model (IBM) was considered. The present work focuses on the fibrous medium only and aims to derive meso and macroscopic models from the background IBM using techniques of kinetic theory. Indeed, the computational cost of an IBM scales polynomially with the number of agents, which makes them practically untractable for large systems. Continuum models allow to break this curse of scaling but they suppose that a suitable coarse-graining procedure which averages out the fine-scale structure has been applied to the IBM. In order to capture the correct effects of the finescale dynamics on the large-scale structures, it is of paramount importance to perform this coarse-graining as rigorously as possible. This is the aim of the present work.
The derivation of a continuum model from the fiber dynamics is done in two steps. We first derive a kinetic model from the underlying IBM and secondly, we perform a diffusion approximation of the latter to obtain the continuum model. The kinetic model provides a statistical mechanics description of the underlying IBM by investigating how the probability distribution of fibers in position and orientation space evolves in time.
Here, we will show that the mere distribution of fibers is not sufficient to close the system and that the cross-link probability distribution needs to be introduced. The cross-links provide correlations between the fibers and consequently their distribution can be viewed as similar to the two-particle fiber distribution. We will formally show that the knowledge of the one-and two-particle distributions is enough to provide a valid kinetic description of the system. Of course, this fact needs to be confirmed by numerical simulations and mathematical proofs. But if it proves correct, this model provides a unique example, to our knowledge, of a kinetic model which is closed at the level of the two-particle distribution function. Indeed, the question whether or not kinetic descriptions must include higher order distribution functions has been actively discussed in the recent years [10,11,25,26] . We also note that the introduction of the cross-link distribution functions provides an economic and efficient way of statistically tracking the fiber network topology. This methodology could prove interesting for other situations of dynamically evolving networks.
The second step consists of a diffusion approximation of the previously derived kinetic model. It starts with changing the time and space units to macroscopic ones. The macroscopic space unit is large compared to the typical spatial scale of the fibers, e.g. their length and the macroscopic time unit is large to the typical time scale of the fibers, e.g. the time needed for two linked fibers to align with each other. A diffusive rescaling relates the time and space rescaling in such a way that the ratio of the microscopic to macroscopic time units is the square of that of the spatial units. This choice is made necessary by the absence of any polarization in the medium which makes diffusive behavior dominate. A key assumption that we make here is to assume that the linking/unlinking frequencies are very large: the typical linking/unlinking time measured in the macroscopic time unit scales like the square of the typical fiber alignment time (also measured in macroscopic unit), which is very small. This allows us to deduce an algebraic relation between the cross-link distribution function and the fiber distribution function, and to realize a closure of the kinetic equation at the level of the fiber distribution function alone. This assumption is questionable given the biological applications we have in mind, but it provides a first step towards a more complete theory involving finite linking/unlinking times.
From these assumptions, we derive a singular perturbation problem for the fiber kinetic distribution function that has the form of a classical diffusion approximation problem [4,16,31], whose leading order collision operator comes from the nematic alignment of the fibers due to the alignment torque at the cross-links. This operator has equilibria in the form of generalized von Mises distributions of the fiber directions. The von Mises distribution extends Gaussian distributions to probabilities defined on the unit circle. It is peaked around a mean fiber direction angle θ 0 . The continuum model describes how the local fiber density ρ and the local fiber direction θ 0 vary as functions of position x and time t. To obtain these evolution equations, we must integrate the kinetic equation against suitably chosen collision invariants. This operation cancels the singularly perturbed term. Here, the difficulty it that there exists only one such collision invariant in the classical sense, which allows us to find an equation for the density ρ only. To find an equation for the mean fiber direction θ 0 , we use the recently developed theory of Generalized Collision Invariants (GCI) [14,15,17,19]. The resulting system is a nonlinear coupled system of diffusion equations for ρ and θ 0 . In the case of a homogeneous fiber distribution, when the density is uniform in space and constant in time, we show that the resulting nonlinear diffusion model for θ 0 is parabolic. In future work, it will be shown that this system is well-posed. Numerical simulations will demonstrate that the continuum model provides a consistent approximation of the underlying IBM for the fiber dynamics. Numerous macroscopic models for fibrous media have been previously considered in the literature but very few of them have been derived from an underlying IBM. Most of them are heuristically derived from continuum theories such as mechano-chemical principles [2,33], thermodynamics [22], or viscous fluid mechanics [23].
The outline of this paper is as follows. In Section 2, we start with the description of the IBM. Section 3 is devoted to the derivation of the kinetic model. The scaling assumptions and the scaled kinetic equations are derived in Section 4. In Section 5, we perform the large scale limit of the so-obtained equations. Finally, Section 6 is devoted to the analysis of the model in the case of a homogeneous fiber density. Conclusions and perspectives are drawn in Section 7. Some technical computations are detailed in Appendices.

Individual Based Model for fibers interacting through alignment interactions
We intend to model a medium consisting of interconnected fibers. To simplify the geometric description of fibers, we decompose them into fiber elements of uniform fixed length and consider that a fiber consists of several connected fiber elements. The link between two connected fibers can be positionned at any point along the fibers (not only the extremities) and a given fiber can be connected to any number of other fibers, thereby allowing to model the branching off of a fiber into several branches. The links are not permanent. The topology of the fiber network is constantly remodelled through link creation/deletion processes. To model fiber resistance to bending, we suppose that pairs of linked fibrs are subject to a torque that tends to align the two fibers with respect to each other. Finally, the fibers are subject to random positional and orientational noises to model the movements of the tissue and to positional and orientational potential forces to model the action of external elements. In the case of a fibrous tissue, these external elements may consist of cells or other tissues.
In this paper, we restrict ourselves to a two-dimensional model. We consider a set of N fiber elements modelled as small line segments of uniform and fixed length L, described by their center X i ∈ R 2 and their angle θ i with respect to a fixed reference direction. As the fiber elements are assumed apolar, θ i is an angle of lines, i.e. θ i ∈ [− π 2 , π 2 ) modulo π. We define energies related to each of the phenomena described above namely an energy for the maintenance of the links W links , an energy for the alignment torque W align , an energy for the action of the external elements W ext , an energy for the noise contribution W noise and a total energy made of the sum of all these energies: All these energies are functions of the N fiber positions (X i ) N i=1 and orientations (θ i ) N i=1 . Note that W noise is rather an entropy than an energy, so that W tot is indeed the total free energy of the system. Fiber motion and rotation during a time interval between two fiber linking-unlinking events is supposed to occur in the steepest descent direction to this free energy, namely according to: To define the expression of W links , we consider a time at which no linking/unlinking process occurs. Then, the set of links is well-defined and supposed to have K elements. Let k ∈ {1, . . . , K} be a given link and denote by (i(k), j(k)) the pair of indices corresponding to the two fibers connected by this link. To make the labeling of the pair unique, we assume without loss of generality that the first element of the linked pair is always the one with lowest index, i.e. i(k) < j(k). The link is supposed to connect two points X k i(k) and X k j(k) on fibers i(k) and j(k) respectively. These points are determined by the algebraic distances ℓ k i(k) and ℓ k j(k) to the centers X i(k) and X j(k) of the two fibers respectively; We thus have the relation: where ℓ k i(k) , ℓ k j(k) ∈ [−L/2, L/2] and where, for any fiber i, we let ω i = (cos θ i , sin θ i ) be the unit vector in the direction of the fiber. All along the link lifetime, the link places a spring-like restoring force that attracts X i(k) back to X j(k) (and vice-versa) as soon as their are displaced one with respect to each other. This restoring force gives rise to a potential energy where κ is the intensity of the restoring force. Obviously, the larger κ, the better the maintainance of the link is ensured. The potential W links is then assumed to be the sum of all the linked fiber spring forces: We stress the fact that the quantities ℓ k i(k) and ℓ k j(k) remain constant throughout the link lifetime. They are determined at the time of the creation of the link (see below and  l ij and l ji refer tol(X i , θ i , X j , θ j ) and ℓ(X j , θ j , X i , θ i ) (2.12). A. Situation at linking time. B. Restoring potential V ij (2.4) after motion of the fibers.
The external potential W ext associated with the external forces is supposed to be the sum of potential forces U(X i , θ i ) acting on each of the N fibers: (2.6) Here, U(x, θ) is a given, possibly time-dependent smooth function. In the case where the system describes the collagen fibers in a tissue, U aims to model the presence of cells or other organs. Linked fibers are subjected to an alignment force at their junction to model fiber resistance to bending. This force tends to align linked fibers i(k) and j(k) and derives from the potential b(θ i(k) , θ j(k) ) which reads: where α plays the role of a flexural modulus and β is a modeling parameter. The binary alignment potential only depends on the angles θ 1 and θ 2 , and the total alignment energy W align is supposed to be the sum of all the binary alignment interactions: We include random positional and orientational motion of the fiber elements which, in the context of tissue dynamics, originate from the random movements of the subject. With this aim, we introduce an entropy term: wheref is a 'regularized density': Here, ξ N and η N are regularization functions which allow to define the logarithm off and have the following properties: and Supp stands for the support of a function. Here, R N and M N are chosen such that √ NR N and NM N → ∞ as N → ∞. The mean interparticle distance in x and θ are respectively of order 1 √ N and 1 N . This condition is equivalent to 1 √ N R N → 0 and 1 N M N → 0, which means that as N → ∞, the number of particles inside the support of a regularizing kernel tends to infinity. This way of modeling the influence of the noise is customary in polymer dynamics [7].In the next section, we show that such an entropy term gives rise to diffusion terms at the level of the mean-field kinetic model. By inserting (2.5), (2.6), (2.8) and (2.9) into (2.2), (2.3), we find the fiber equation of motion, during any time interval between two linking/unlinking events: which we can write: When two fibers i and j intersect each other, because of the continuity of their motion, they are going to intersect each other during a time interval [t * , t * ]. We assume that, during this time span, the linking probability follows a Poisson process of parameter ν f , i.e. the probability that a link is formed during the interval [t * , t] with t < t * is 1 − e −ν f (t−t * ) .
Only one link can be formed between the two fibers of the same fiber pair. Supposing that a link, indexed by k is formed between the fibers i and j (such that i = i(k) and j = j(k) if i < j) at a time t k ∈ [t * , t * ], we denote by X k the attachment site of the link. The distancel(X i(k) , θ i(k) , X j(k) , θ j(k) ) between the center X i(k) of fiber i(k) to the k-th link attachment site X k with fiber j(k) (see Figure 1.B) can be directly computed by: ) are the coordinates of the center of fiber i(k). For X = (x, y) and ω = (α, β), we denote by X × ω = xβ − yα. Then,l(X i(k) , θ i(k) , X j(k) , θ j(k) ) can be written:l where again, ω(θ) = (cos θ, sin θ) is the directional vector associated to angle θ. The fact that the two fibers are intersecting each other at time t k is written: where L is the fiber length and where all positions and angles are evaluated at time t k . The quantitiesl( ) at the time t k of the formation of the link set the positions of the attachment sites X k i(k) and X k j(k) of the link on fibers i and j. Therefore, ℓ k i(k) and ℓ k j(k) remain constant throughout the link lifetime and equal to their value at the time t k . So, we have throughout the lifetime of the link. We also assume that existing links can disappear according to a Poisson random process of parameter ν d , i.e. the probability that the link disappears in the time interval The next section is devoted to the asymptotic limit N, K → ∞ of this model.

Derivation of a kinetic model
Here, the derivation of a kinetic model from the Individual Based Model of section 2 is performed. The empirical measure f N (x, θ, t) of the fibers is introduced: where δ (X i (t),θ i (t)) (x, θ) denotes the Dirac delta located at (X i (t), θ i (t). It gives the probability to find a fiber at point x and orientational angle θ at time t. The empirical measure g K (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 , t) of the fiber links is given by: , with a similar definition of the Dirac deltas. It gives the probability of finding a link with associated lengths within a volume dℓ 1 dℓ 2 about ℓ 1 and ℓ 2 , this link connecting a fiber located within a volume dx 1 dθ 1 π about (x 1 , θ 1 ) with a fiber located within a volume dx 2 where f and g satisfy equations given in the following theorem: and S(g) is given by: where δl(ℓ 1 ) denotes the Dirac delta atl, i.e. the distribution acting on test functions This kinetic model consists of two evolution equations. The first one (Eq. (3.1)) is an equation for the individual fibers and describes the evolution of the one-particle distribution function f . Eq. (3.2) is an equation for the links between fiber pairs. The distribution function g describes the correlations between fiber pairs brought by the presence of links. It can be viewed as a kind of two-particle fiber distribution function. This model is, to our knowledge, a unique explicit example of a kinetic model written in terms of the one and two particle ditributions and closed at this level. Also, the distribution function g can be seen as a way of describing the random graph of the fiber links, namely the graph where the nodes are the fibers and the edges are the links. This statistical description of a random graph could be useful to describe other kinds of random networks, notably in social sciences. As the links are tightly tied to the fibers, they are convected by them and follow their motion. Simultaneously, they constrain the linked fibers to move together, so they directly influence their motion. The action of the links on the individual fiber motion is contained in the third and sixth force terms F 1 and F 2 of Eq. (3.1) and are the kinetic counterparts of (2.4). The second and fith terms describe transport in physical and orientational spaces due to the external potential and are the kinetic counterparts of (2.6). The fourth and seventh terms are diffusion terms of amplitude λd and µd respectively. They represent the random motion of the fibers and originate from the interactions described by Eq. (2.9). The individual motion of the fibers is thus related to the motion of its linked neighbors. The left-hand side of Equation (3.2) describes the evolution of the links between fibers. Indeed, it is composed of the convective terms generated by the external potential and by the diffusion terms. The forces induced by the restoring potential generated by the links again gives rise to the non local terms F 1 and the first term of F 2 . The kinetic counterpart of the alignment force between linked fibers (see Eq. (2.8)) is encompassed in the second term of the force F 2 and only acts on the orientation of the fibers. The right hand side S(g) of equation (3.2) describes the Poisson processes of linking/unlinking at frequencies ν f and ν d , respectively. The first term describes the formation of the link and the Dirac deltas indicate that, at the link creation time, the link lengths ℓ 1 and ℓ 2 are set by the geometric configuration of the fibers at the attachment time. Also, because ℓ 1 and ℓ 2 are restricted to lie in the interval [−L/2, L/2], we see that the link creation term is non-zero only when two fiber elements are intersecting each other. The second term just describes a decay of the link distribution at the rate set by the Poisson process, i.e. ν d .
The formal proof of this result is inspired from Ref. [32], and the detailed computations can be found in appendix A. The rigorous proof of this result is an open question and is left for future work.

Dimensionless Equations
We express the problem in dimensionless variables. For this purpose, let t 0 be the unit of time and x 0 , f 0 = 1 the units of space, distribution function and energy. The scaling of f (x, θ) and g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) comes from the fact that they are probability distribution functions on a 2D domain. The following dimensionless variables are defined: . and the following dimensionless parameters are introduced: First of all, from the expression of V (see Eq. (2.4)), we get: , one notes that: . Finally, if the space and time scales x 0 , t 0 are chosen such that λ ′ = χ = 1, i.e: the dimensionless equations forf andḡ read (dropping the primes and tildes for the sake of clarity):

Scaled equations
So far, the chosen time and space scales are microscopic ones, and describe the system at the scale of the agent interactions. In order to describe the system at a macroscopic scale, a small parameter ε ≪ 1 is introduced and the space and time units are set tõ The fiber length measured at scale x 0 is supposed to stay of order 1 as ε → 0, i.e. L = O(1). The variables x, t, ℓ and unknowns f and g are then correspondingly changed tox = where U 0 is acting on the space variable only and U 1 is a π-periodic potential acting on fiber orientation angles only. The external potential acting on the space variables is supposed to be one order of magnitude stronger than the one acting on the fiber rotations: The strength of the alignment potential is supposed to be large α = O(ε −1 ), i.e.α = εα withα = O(1), and we choose the exposant β = 1. The intensity of the alignment potential between linked fibers is supposed to be small κ = O(ε), i.e. κ = ε −1 κ withκ = O(1) and the diffusion coefficient and parameter ξ are supposed to stay of order 1: d, ξ = O(1). In order to simplify the analysis of the system, the process of linking/unlinking is supposed to occur at a very fast time scale, i.e.ν f = ε 2 ν f and and consequently, Then we have: Finally, we define X 1 and X 2 such that: The macroscopic fiber linking/unlinking operator S(g) is similar to the one defined Eq. (3.5). Indeed, from Eq. (2.12):l(x 1 , θ 1 , x 2 , θ 2 ) = ε −1/2l (x 1 , θ 1 ,x 2 , θ 2 ) and thus: Altogether, the macroscopic version of Eqs. (4.1)-(4.2) reads (dropping the tildes for the sake of clarity): From now on, we note f ε =f and g ε =g. The following proposition holds: Assuming f ε and g ε exist, then, formally, they satisfy:

5)
and and Remark 4.1. In the proof of proposition 4.1, we will show that (4.11) (4.12) The proof of this proposition is given in section 4.3. From these equations, one notes that the hypothesis of dominant creation/deletion of links makes the reaction forces F 1 and F link of order O(ε 3 ). In this case, the process of linking/unlinking is so fast that the constraint is satisfied at all times. Moreover, under this assumption, the first contribution of the alignment force acting on a fiber is the sum of elementary alignment forces generated by its intersecting fibers, weighted by ν f ν d . One also notes that the alignment force F al is local in space.
Under these scaling assumptions, the leading order of the left-hand side of Eq. (4.5) takes the form of a collision operator of kinetic theory. It acts on the orientation vector θ only and it expresses that the alignment potential (2.8) is counter-balanced by the diffusion term which tends to spread the particles isotropically on the sphere. The other terms act at lower order ε.
As the large scale limit involves an expansion of the solution around a local equilibrium, the study of the local equilibria of the collision operator are of key importance. Therefore, section 5 will be dedicated to the study of the properties of the left-hand side of (4.5).

Large scale limit
In this section, the limit ε → 0 of the solution f ε to (4.22) is explored. For this purpose, Eq. (4.22) is rewritten where the collision operator Q(f ε ) is defined by and where we recall that C 1 and G[f ] are defined by (4.8) and (4.9) respectively. The operator Q is a non linear operator on f which acts on θ only and leaves x and t as parameters. For each function Φ(θ), we define M Φ (θ) by: where Z is a normalization factor such that Z = π 2 − π 2 e −ξΦ(θ)/d dθ π . Thus, M Φ (θ) is a probability distribution of θ. Such functions are called generalized Von Mises distributions (the Von Mises distribution being the case of Φ(θ) = − cos θ). The next section is devoted to the analysis of the properties of Q(f ) and follows closely Ref. [13].

Equilibria
In this section, the equilibria of the operator Q are studied, and the following proposition is proven: Proposition 5.1. Here, we restrict ourselves to functions of θ only. (i) The operator Q can be written: ) . This proposition shows that the equilibria of operator Q are generalized Von Mises distributions of θ, weighted by the particle density.
Proof. To prove (i), one can note that: To prove (ii), note that f = ρM Φ[f ] is solution of (5.5). Conversely, suppose that f is such that We define the sets H f and V f by: The norms · H f , · V f on H f and V f are then defined such that: and For f ∈ V f using Green's formula, we get: , with ρ > 0, which ends the proof. Now, the following lemma is proven: For any function f (θ), the potential function Φ[f ](θ) of Eq. (5.3) can be written: where C 1 is given by (4.9), C = or equivalently by: Remark that the second condition is equivalent to saying that and this defines θ f uniquely modulo π.
Thanks to Eq. (5.11), Hypothesis 5.1 amounts to supposing that the ratio ν f ν d is inversely proportional to the fiber density.
Since there is no obvious conservation relation other than the conservation of the local fiber density, the only collision invariants in this model are the constants. The integration of equation (4.5) against these invariants does not allow us to find the evolution equation for the mean orientation. In order to obtain an equation on θ 0 , inspired from Ref. [17], the concept of Generalized Collision Invariants (GCI), i.e. of collision invariants when acting on a restricted subset of functions f , is introduced.

Collision invariant
A collision invariant is a function Ψ such that for all function f of θ, Q(f )Ψdθ = 0. However, due to the lack of momentum conservation, the only collision invariants are the constants. This is not enough to determine both ρ and θ 0 . To this aim, following Refs. [19] and [17], we introduce the notion of GCI. For any θ 0 ∈ [− π 2 π 2 ), we define L θ 0 as the following linear operator: Definition 5.4. For a given θ 0 ∈ [− π 2 , π 2 ) a GCI associated to θ 0 is a function Ψ such that: The set of the GCI associated to a given θ 0 ∈ [− π 2 , π 2 ) is a linear space denoted by G θ 0 . Lemma 5.5. Ψ ∈ G θ 0 if and only if ∃β ∈ R such that: where L * θ 0 is the L 2 formal adjoint of L θ 0 , i.e.
Proof. By (5.7), the condition θ f = θ 0 mod(π) is equivalent to the linear constraint: By a classical duality argument [17], we deduce that Ψ ∈ G θ 0 if and only if: Note that now, there are no more constraints on f . Therefore, we can eliminate f and get (5.14).
Equation for θ 0 We multiply Eq. (5.1) by the GCI Ψ θ f ε associated with the direction θ f ε of f ε , namely Ψ θ f ε = g(θ − θ f ε ) where g is the function defined in Prop. 5.6. We integrate with respect to θ and first note that: by (5.13). Since f ε → ρM θ 0 , we have θ f ε → θ 0 and Ψ θ f ε → Ψ θ 0 . Therefore, in the limit ε → 0, we get: (5.24) For simplicity, we denote M θ 0 = M. We have: Using the continuity equation (5.20), we have: So: Therefore, Eq. (5.24) reads: where: We now turn to the development of each term of Eq. (5.25). We have: Then, and thus, X 1 can be written: From integration by parts, the following relations can be written: Therefore, we have: Since X 2 is the integral of a π-periodic function over a period, we can write Now, by construction, (see prop 5.6), Ψ(θ 0 − π 2 ) = Ψ(θ 0 ) = Ψ(θ 0 + π 2 ) = 0. So, integrating by parts, we have Now, by construction again (see (5.18)), we have Using again the π-periodicity of U 1 , we obtain: Now, let us turn to X 3 . The details of this computation are postponed to appendix B. We find: where, using (5.8), We note that αL 4 γ 48r γ 1 = 1 2r α 3 . Finally, let us explicit the last term X 4 . A direct computation gives: Then, we deduce that By symmetry, we have: Therefore, with (5.31), we get: Collecting (5.32) to (5.36) and inserting them into (5.25) leads to (5.21).

Case of a homogeneous fiber distribution: stationary solutions
In this section, we study the stationary solutions of (5.20)-(5.21) in the case of a spatially homogeneous fiber distribution and consequently no external spatial potential U 0 = 0. We make the following assumption: Hypothesis 6.1. The fiber spatial distribution is supposed to be homogeneous, i.e. there exists a constant ρ 0 > 0 such that ρ(x, t) = ρ 0 for all (x, t) ∈ R 2 ×[0, ∞). We also suppose that there are no external spatial forces, i.e. U 0 = 0.
We first note that in the absence of external forces, a uniform and constant density ρ 0 is a solution of Eq. (5.20). Now, we are interested in the stationary solutions for the fiber orientation equation (5.21). Noting that the terms involving the spatial derivatives of ρ, we find that such stationary solutions satisfy the following equation: In this equation, the coefficients r, α 1 , α 2 and α 3 are constants thanks to (5.8). Moreover, using (5.22), they can be written as functions of d, L 2 and r as follows: We now show that (6.1) is an elliptic equation. We first introduce some definitions.
Given a function f(x, E) smooth in its arguments x ∈ Ω, E ∈ R × R 2 × S 2 (R), where S 2 (R) is the space of 2 × 2 symmetric matrices with real coefficients, we define the non linear differential operator F : C ∞ (R 2 ) → C ∞ (R 2 ) such that for any x ∈ R 2 and any . The operator F is said to be elliptic at u 1 ∈ C ∞ (R 2 ) (see Ref. [34]) if its linearization DF (u 1 ) is an elliptic, linear differential operator. We state the following proposition: Proposition 6.1. Eq. (6.1) can be put in the form where f(x, D 2 θ 0 ) is the following operator, quasi linear in θ 0 : Here, h(θ 0 ) = ∂ θ U 1 and A(θ 0 ) = (a ij (θ 0 )) i,j=1,2 is a 2 × 2 matrix such that: Moreover, if the following condition is satisfied for all r ∈ R + : where A(r) is given by (6.5), then F (θ) = f(x, D 2 θ) is elliptic at θ 1 for all θ 1 ∈ C 2 (R 2 ).

Conclusion
In this paper, we have formally derived a macroscopic model for temporarily linked fibers interacting through alignment at the links. We have shown that the corresponding kinetic model involves two distribution functions: the fiber distribution function and the crosslink distribution function. The latter can be seen as a joint two-particle fiber distribution function. This model provides a unique explicit example of a kinetic model closed at the level of the two particle distribution function. We then considered the regime of a fast fiber linking/unlinking process, where the link distribution function can be expressed simply in terms of the fiber distribution function. We studied the diffusive limit of the resulting equation and obtained a system of two coupled nonlinear diffusion equations for the fiber density and mean orientation. In the homogeneous fiber density case, we showed that the resulting quasilinear problem is elliptic. Future works will deeper investigate the mathematical properties of the models, such as rigorously proving the mean-field kinetic limit of the particle system or proving existence and uniqueness of smooth solutions for the macroscopic diffusion system. Numerical simulations will be performed to validate the macroscopic model by comparison with the individual based model. Further perspectives are the removal of the fast fiber linking/unlinking hypothesis, in order to understand how a finite lifetime of the cross-links affects the macroscopic dynamics.
A Proof of Theorem 3.1

A.1 Evolution equation for the fibers
For all observable functions Φ(x, θ), we define: Similarly, for all two-particle observable functions Ψ(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ), we define: where integrals over x are carried over R 2 , in θ over (− π 2 , π 2 ) and in ℓ over (− L 2 , L 2 ). We recall the notations Using (2.10) and (2.11), we obtain: We get, using the definition of a distributional derivative: Now, exchanging the sums in i and k in the previous equation, one obtains: From the symmetry of V (see Eq. (2.4)), the following expressions hold: and from the symmetry of b, we have: leading to: or again: Therefore, we obtain: Finally, we get: where, for a distribution T acting on functions of (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ), we denote by [[T ]](x 1 , θ 1 ) the distribution which to any function Φ(x 1 , θ 1 ) associates and where 1 is the constant function of the variables (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) equal to 1. In the formal limit N → ∞, K N → ξ and given the assumptions on the regularizing sequences θ f and we obtain: where,

A.2 Evolution equation for the fiber links
Following the same principle as for f N and given that the links are maintained over time, i.e.
, one can write: where E k corresponds to the k-th line of (A.3). For the sake of simplicity, the computation of E 1 only is developed here. The computation of the other ones are similar and omitted. From Eqs. (2.2), (2.3), one obtains: where we write V = V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ). Now, exchanging the sums in k and k ′ and using the symmetry of V , one obtains: Because there is no restriction on the number of links per fiber, the sums over k cannot be simplified in this case. In order to express the third and fourth terms, the number C k ′ i (resp. C k ′ j ) of fibers linked to fiber i(k ′ ) (resp. j(k ′ )) is introduced: where Card denote the cardinal of a set. Then, as K → ∞, the following expression holds for any chosen fiber k ′ : is the conditional probability of finding a link conditioned on the fact that one of the fibers of this link has the same location and orientation as i(k ′ ). Then, as N → ∞, K → ∞ such that K N → ξ > 0 , C k ′ i is the mean number of links per fiber. The mean number of links in the volume π dℓ 1 dℓ 2 and the mean number of fibers in dX i(k ′ ) dθ i(k ′ ) is Nf (X i(k ′ ) , θ i(k ′ ) ). Thus: .