A Novel Modeling and Simulation Approach for the Hindered Mobility of Charged Particles in Biological Hydrogels

This article presents a novel computational model to study the selective filtering of biological hydrogels due to the surface charge and size of diffusing particles. It is the first model that includes the random 3D fiber orientation and connectivity of the biopolymer network and that accounts for elastic deformations of the fibers by means of beam theory. As a key component of the model, novel formulations are proposed both for the electrostatic and repulsive steric interactions between a spherical particle and a beam. In addition to providing a thorough validation of the model, the presented computational studies yield new insights into the underlying mechanisms of hindered particle mobility, especially regarding the influence of the aforementioned aspects that are unique to this model. It is found that the precise distribution of fiber and thus charge agglomerations in the network have a crucial influence on the mobility of oppositely charged particles and gives rise to distinct motion patterns. Considering the high practical significance for instance with respect to targeted drug release or infection defense, the provided proof of concept motivates further advances of the model toward a truly predictive computational tool that allows a case- and patient-specific assessment for real (biological) systems.


Introduction
The remarkable ability of hydrogel forming biopolymer networks to control the mobility of different kinds of diffusing molecules and particles individually is of crucial importance for numerous functions of the human body. See [1] for a recent review of this topic. On the one hand, this selective permeability gives rise to the protection of the organism against pathogens such as viruses that are effectively hindered from invading and traversing the organism. On the other hand, it ensures the effective transport of a broad variety of substances that are useful and important for the organism. Examples of such biological hydrogels include mucus, the extracellular matrix (ECM), intracellular biopolymer networks (comprising, e.g., neurofilaments or actin), the vitreous humor and the matrix of biofilms, and can thus be found throughout the entire human body and in countless other places in nature. The high practical relevance of this fundamental topic extends over multiple fields ranging from medical diagnosis to the therapy of body malfunctions and targeted drug release. This creates yet unimagined possibilities in technical applications such as filters used in chemical, mechanical or medical engineering.
A large number of experimental studies have investigated the origin of this selective filtering. Meanwhile, there is strong evidence that -besides the most obvious mechanism to filter by size -the surface properties of the particles also plays an important role (see e.g., [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]). This new paradigm thus suggests filtering not (only) by size but also through a combination of other particle properties such as charge, hydrophobicity or binding affinity and can thus be referred to as interaction filtering. Despite considerable scientific effort in this field, many aspects concerning the underlying mechanisms and specific conditions remain unknown. To a large extent, this can be explained by the following three factors. First, there is a great complexity in the many different biological systems both in number and diversity of components, e.g., with respect to their molecular architecture. Second, there are big challenges and limitations with respect to experimental preparation and measurement techniques when it comes to the required high spatial and temporal resolution, especially over considerable time spans of several seconds. And third, there is a fundamental lack of in-depth understanding of physical and chemical interactions on the molecular scale. This lack of detailed, fundamental microscopic understanding prevents any reliable prediction of the diffusive mobility of molecules and particles in other than the few particle-hydrogel systems already studied in vitro. At this point, the development of accurate and efficient computational models, which are capable of resolving small scales and covering large spans in space and time, is expected to substantially contribute to scientific progress in this field.
Compared to the large body of literature reporting on experimental work in this area, relatively few computational studies have been published so far. Most of the early and also some of the recent work focused solely on the excluded volume effect of the fibers, i.e., the repulsive steric interactions with the network fibers that hindered the free diffusion of a (hard) spherical particle (see e.g., [17][18][19][20][21][22]). In his Monte Carlo simulations, Saxton [23] was the first to include and study other than steric interactions in the form of a binding model. Since the recognition of the dominant role of electrostatic and possibly other types of molecular interactions as outlined above, several computational studies have included electrostatic effects and confirmed the trends observed in experiments and shed light on the underlying mechanisms [13,[24][25][26][27][28]. Particularly the recent works published by Hansing et al. [27] and Hansing and Netz [28] were very successful in analyzing the influences of particle size, fiber volume fraction, particle charge, and the comparison of oppositely vs. similarly charged particles and networks. They also found good agreement of the simulation results with several sets of experimental data, which confirms the validity of such modeling approach.
The computational model proposed in this work aims to improve especially the modeling of the fiber network, which has been modeled in a very simplified manner in all previous studies. Either a square array of straight and parallel rigid fibers oriented along one spatial dimension [25] or a cubic lattice consisting of either linearly aligned hard spheres [26] or consisting of straight rigid fibers [13,27] has been assumed. Zhou and Chen [24] likewise applied a cubic lattice consisting of beads placed at the vertices and connected by linear spring elements. Hansing and Netz [28] were the first to break the strong geometrical order in all of these models, they only allowed, however, for a variation of the spacing of the still straight and infinitely long rigid fibers with a mutually orthogonal orientation. Yet they found a significant and fundamental difference in the trapping mechanism of ordered and disordered fiber lattices, which can be attributed to locally denser regions of the network in the case of attractive particle-network interactions. This is a strong motivation to work toward a more realistic modeling of the fiber network as a crucial part of biological hydrogels as will be outlined in the following paragraphs.
In our approach, the individual fibers will be modeled by the geometrically exact 3D beam theory, thus allowing for arbitrarily curved initial shapes of the fibers and possibly large deformations, including all six modes of axial strain, (2x) shear, torsion and (2x) bending. In addition, the initial spatial distribution, orientation and interconnection of fibers will be modeled as a random 3D Voronoi network, thus mimicking several important geometrical features of real biopolymer networks such as a random, spatially variable mesh size distribution, a random fiber length distribution, random mutual orientations of the fibers, and arbitrary connectivity between the fibers. Several of these attributes are expected to play an important role when it comes to both purely steric interactions and the combination with electrostatic interactions and will be studied in Section 3. A similar approach to network generation based on Voronoi tessellation has previously been applied in a number of publications, e.g., to study cell-cell communication in a 2D network of linearly elastic springs [29].
Applying such a sophisticated model to the individual fibers and the network they constitute comes at the cost of an increased complexity and size of the system of equations to be solved in the simulations. However, the same modeling strategy based on geometrically exact beam finite elements describing the biopolymer fibers has previously been applied to model the Brownian dynamics of individual semiflexible filaments [30,31] and has been proven to be accurate and also efficient enough to study large-scale problems such as the process of self-assembly of (different morphologies of) transiently cross-linked biopolymer networks [32,33] as well as their (high-and low-frequency) rheology [34]. On the one hand, this confirms that the novel modeling approach is suited to study the hindered mobility of particles in hydrogels and on the other hand, this already outlines the long-term opportunities. Using the self-assembled network configurations, e.g., for an actin bundle network, can readily replace the Voronoi-type network applied as a first step in the present study. Moreover, the dynamics of the network, including the reorganization of the transient cross-links could be included and would be highly interesting since there is experimental evidence that particles larger than the mesh size can still diffuse through the network by locally breaking inter-fiber links [6]. This has also been confirmed to be an effective transport mechanism in a first computational model [35]. Another example for the dynamics of hydrogel networks is given by the self-renewal process of a mucus layer [1,36].
In general, the mathematical description of the problem and therefore also the numerical methods required to solve it as well as the code framework substantially differs from the ones used by the previous computational studies listed above. This is an inevitable consequence of the aforementioned modeling of individual biopolymer fibers as elastically deformable, geometrically exact beams. In the present model, a set of well-established numerical formulations for beams [37,38] is combined with novel beam-sphere interaction models, which are derived as a special case of the more general approach to fiber-fiber interactions recently proposed in [39]. The latter is the first 3D beam-beam interaction model for molecular interactions between curved slender fibers undergoing large deformations and is thus an important prerequisite for the beam-sphere interaction formulations to be applied in this article.
The remainder of this article is structured as follows. Section 2 presents all aspects of the computational model and the required numerical methods. After the presentation and discussion of the results in Section 3, this article will be concluded with a summary of the findings and an outlook to promising aspects of future research in Section 4.

Computational model and methods
Given the broad variety of biological hydrogels mentioned above, we chose the ECM gel(s) used in the comprehensive experimental studies of particle mobility by Lieleg et al. [7] and Arends et al. [12] to serve as the main reference for the specific setup and parametrization of the versatile computational model. For these gels, an additional, subsequent characterization in terms of their biophysical properties has been conducted [40], which shall prove useful in the following. To reduce the complexity of this multi-component biological system, the computational model considers only the following three key components, which are thought to mimic the crucial influences of the studied problem: a fiber network, a diffusing particle, and their mutual interactions. All three model components will be considered individually in the following and selected further aspects of the simulation setup such as boundary conditions and post-processing of the results will be discussed. The software package used for the simulations in this article is the parallel, multi-physics, in-house research code BACI [41].

(a) Overall approach
Following up on the concept of modeling biopolymer fibers as elastically deformable, geometrically exact beams as described above, the overall approach follows the one commonly used in nonlinear finite element frameworks for structural dynamics. In short, this approach consists of the following steps: According to the principle of virtual work, the weak form of the mechanical balance equations is derived and subsequently discretized in space and time. Given a proper set of initial and boundary conditions, an implicit load/time stepping scheme is applied and in every step the solution of the resulting discrete system of nonlinear equations is found iteratively by means of Newton's method. Refer to Ref. [39] for a discussion of selected aspects of the applied algorithms and libraries.
(b) Biopolymer fiber network As outlined above, the initial, stress-free configuration of the fiber network is the result of a 3D Voronoi tessellation of the cubic simulation box (size 10 × 10 × 10 µm) generated via the open source library voro++ [42]. Figure 1(a) shows an example of the resulting network architecture. The main input of this preprocessing step are the randomly chosen locations of a number n VP of so-called Voronoi points 1 . The output are the vertices and edges of a random, irregular, polygonal network that are used to define the position and orientation of the initially straight beam segments as well as their interconnections. For Voronoi-based tessellation, the connectivity number, i.e., the number of fibers branching off at each junction point, is 4, which agrees well with values of 3 to 4 reported for ECM gels [44]. In the future, this parameter could also be adapted by randomly removing or adding fibers to match other hydrogel architectures. At each junction, the beam endpoint positions and rotations are coupled, which corresponds to a model assumption based on the expected mechanical behavior of the chemical binding between the fibers, and could be once again adapted e.g., to hinged connections in a straightforward manner. Also, the binding and unbinding of such connections at given rates could be included by following the approach described e.g., in [34]. To be able to use the simulation box as a representative volume element (RVE) with periodic boundary conditions (see Section (f) for details), the Voronoi tessellation and thus the resulting network geometry is chosen to be periodic in space. Altogether, a simple visual comparison with electron microscopy images of real biological hydrogels in Figure 1 reveals a high similarity both in terms of the random, irregular, polygonal structure as well as its characteristic properties, such as the distribution of pore sizes, fiber segment length and connectivity.
The individual fibers are modeled by (geometrically exact) 3D beam theory, assuming undeformable cross-sections of circular shape and a hyperelastic material law. Specifically, we 5 rspa.royalsocietypublishing.org Proc R Soc A 0000000  apply the Simo-Reissner beam theory [45][46][47], which accounts for 6 deformation modes of axial tension, (2x) shear, torsion and (2x) bending. Further specifications such as the dimensional and constitutive parameters are chosen to mimic collagen I as the key constituent of the targeted class of ECM hydrogels, however, all these parameters can be easily adapted to study their influences or to model other fiber species. In the present study, the cross-section diameter is set to D f = 75 nm, Poisson's ratio is set to ν = 0.3, and Young's modulus is varied from as low as E = 0.1 MPa up to the theoretical limit of rigid fibers to study the influence of the fiber stiffness. Note that experimental measurements for collagen I suggest values in the wide range of E = 1 MPa − 1 GPa [48,49], which are covered in this work as well. Based on the results in Section (c), it will turn out that the influence of the fiber stiffness on the particle mobility in the problem setup considered here is negligible for the realistic range of values for Young's modulus E, and a noticeable difference in results can only be observed below a threshold value of E * = 1 MPa. It can therefore be assumed that the specific values for the fiber material are of minor importance and that the deformation of fibers will only become important if the dynamics of the network reorganization will be included as outlined above, or if much thinner and softer, i.e., much more compliant fibers are considered e.g., in the context of a different kind of hydrogel (e.g., mucin or F-actin) or a dysregulation of fiber stiffness.
In order to characterize the generated fiber networks, Figure 2 shows the resulting values of the mean and standard deviation of the fiber volume fractionV f and the minimum/average/maximum cell diameters 2 as a function of the number of Voronoi points n VP . For this statistical analysis of the random network geometries, a box size of 10 × 10 × 10 µm and a fiber diameter D f = 75 nm as stated above is assumed. A number of n VP = 60 Voronoi points result in a fiber volume fractionV f ≈ 0.4% and a range of cell diameters of approx. 2.7 − 5.2 µm. This turns out to match the typical mesh size of 2 − 3 µm reported in [7] quite well and is thus chosen as the default value for most of the simulations conducted in this work (once again refer to Figure 1 for a visual comparison of model and real hydrogels). The densest network to be considered in this work is given as n VP = 600 and thus results in a high fiber volume fraction ofV f ≈ 1.8% and cell diameters in the range of approx. 1.0 − 2.5 µm. This second value of the network parameter n VP = 600 is motivated by the example of the human amniotic basal membrane [44], which appears to be much denser than the one considered above. One example of the resulting model network is shown in Figure 3. In order to investigate and average out the influence of the specific network geometry, 5 different realizations of the random network generation process will be considered as input for the simulations for each -otherwise identical -set of input parameters in Sections (b) and (c), respectively. In contrast to the spherical particle, which will be discussed next, the fibers are assumed to be athermal because the lengths of the individual fiber segments are generally much smaller than their persistence length p ≈ 39 µm (resulting from the parameters for the 7 rspa.royalsocietypublishing.org Proc R Soc A 0000000 (b) (c) Figure 3. (a) Image of an example network resulting for a second, very high value for the fiber volume fractionV f ≈ 1.8% considered in the simulations (box size is 10 × 10 × 10 µm, rendered using Blender [43]). Comparison of (b) a magnified part of this network with (c) a magnified part of the network with the default fiber volume fractionV f ≈ 0.4% (shown in Figure 1(a)).
bending stiffness given above and the thermal energy at room temperature k B T ≈ 4 × 10 −3 aJ) and thermal undulations will therefore be negligible. However, for another network architecture with longer fiber segments or a different species of thinner and thus more flexible fibers such as mucin or F-actin, this effect might become important and may be included in the novel computational model by means of the formulation proposed in [30]. Finally, the spatial discretization of the fibers makes use of the geometrically exact Hermitian Simo-Reissner beam elements proposed in [37], which are mainly based on the wellestablished element formulation proposed by Crisfield and Jelenić [50,51]. The applied centerline interpolation using cubic Hermite polynomials ensures both high accuracy in terms of spatial approximation and a C 1 -continuous geometry representation. This is particularly important for smooth contact kinematics and smooth interaction force distributions as has been shown in the context of beam-beam interactions in [39,52]. When creating the finite element discretization for each of the random network configurations, a default element length l ele = 1 µm is used for each straight fiber segment and one additional shorter beam element is created for the remainder of the random segment length if needed. Based on the previous experience (from challenging scenarios with large deformations such as presented in [39]) with this kind of beam finite element featuring fourth order spatial convergence, this is considered to be a sufficiently fine spatial resolution for the magnitude of deformations observed throughout this computational study, such that the influence of the spatial approximation error on the results is expected to be negligible. To conclude, this strategy typically leads to a number of fiber segments n f ≈ 800, a number of beam elements n ele ≈ 1600, and a number of nodes n node ≈ 4000 for n VP = 60. Likewise, for the densest network with n VP = 600, we obtain n f ≈ 6200, n ele ≈ 6200, and n node ≈ 18700.
(c) Spherical particle The particle is modeled as a rigid sphere and is therefore uniquely described by its midpoint position rp(t) ∈ R 3 as a function of the time t as well as its diameter Dp = 1 µm that is chosen to comply with the experimental studies presented in [8,12]. As in several of the aforementioned previous computational studies (e.g., [27,28]), the Brownian dynamics of the particle is modeled by the Langevin equation (see e.g., [53]), including the stochastic thermal as well as viscous drag forces that are repeated here for the reader's convenience. The velocity-proportional drag force f s,visc implicitly models a quiescent surrounding fluid and makes use of the friction coefficient of a sphere γ = 3πηDp according to Stokes: Here,ṙp denotes the particle velocity, and the fluid viscosity is chosen as η = 1 mPa s, which corresponds to the viscosity of water at room temperature and has been found to be in good agreement with experimental tracer studies of freely diffusing particles [12]. Following the consistent modeling of the Brownian dynamics of slender biopolymers within the space-and time-discrete theoretical framework of the nonlinear finite element method according to [30], the stochastic thermal forces acting on the sphere are given as: Here, the thermal energy is set to k B T = 4.1 × 10 −3 aJ corresponding to room temperature and the last term describes the space-time white noise resulting from a two-dimensional Wiener process W(x, t). Both contributions from viscous and stochastic forces on the spherical particle are added to the total virtual work of the system, which other than that, includes the contributions from internal elastic forces of the fibers and the contributions from the particle-fiber interactions that will be discussed next. Afterward, the temporal discretization and time stepping scheme will be presented in Section (e).

(d) Interactions between the particle and the fiber network
In view of the central research questions described above, the interactions between the spherical particle and the fiber network are the key component of this computational model. In accordance with the sophisticated modeling of elastic fibers, these interactions are modeled at the level of individual particle-fiber pairs and are evaluated as a fully resolved, resulting line force distribution on the fiber. At this point, it becomes clear that the sophisticated fiber model, including potentially large 3D deformations, which is a unique feature of this novel computational model, carries over to interaction modeling in the form of additional challenges to accurately and efficiently describe the sphere-fiber interactions for arbitrarily deformed fiber configurations and mutual orientations. Since the conclusion of previous experimental as well as computational studies is that the combination of repulsive steric and (attractive) electrostatic effects is the main reason for the effective selective filtering of biological hydrogels (see e.g., [8,12,28]), both interaction types are accounted for in the present model.
To the best of the authors' knowledge, the problem of a beam interacting with a rigid sphere via both repulsive steric and (attractive) electrostatic effects has not been considered in the literature before. The key idea of our modeling approaches is to consider the sphere-beam interaction as a special, simpler case of the beam-beam interaction. Starting from the formulation for molecular interactions proposed in [39] and the contact formulation proposed in [54], the following procedure basically replaces one of the beams of an interaction pair by a rigid sphere. The resulting formulations, derived for both contact and electrostatic interactions between a rigid sphere and a beam, are summarized below. Note that hydrophobic interactions are not included here due to the lack of a proper fundamental understanding and modeling approach. Assuming that this effect can be described by an effective interaction potential, it could however be incorporated in a very similar manner as the electrostatic interactions and would be an interesting future extension of the present model.

(i) Electrostatic interactions
Generally, the two-body interaction potential of a beam-sphere pair Π ia shall be described as: where a section-sphere interaction potential lawπ section-sphere has been introduced. Such a reduced interaction law is an analytical description of the effective interaction of one cross-section of the beam with the rigid sphere and will be further specified below. In general,π section-sphere will depend on the relative distance vector r b-s of the section midpoint r b (s) and sphere midpoint rs as well as the relative orientation expressed by the relative rotation vector ψ b-s . In order to arrive at the total two-body interaction potential Π ia , the section-sphere interaction potentialπ section-sphere is then numerically integrated along the arbitrarily deformed centerline curve of the beam. Here, s ∈ [0, l] denotes the arc-length parameter, which is defined in the stressfree, initial configuration of the beam's centerline curve. This general approach is in close analogy to the concept of section-section interaction potential (SSIP) laws introduced in [39].
Since it is the most relevant specific example for this study, we now consider the case of Coulomb interaction. Again in close analogy to the careful choice of the SSIP law for this kind of interaction in [39], the following section-sphere interaction potential lawπ section-sphere is obtained: Here, λ b denotes the line charge density of the beam, Qs denotes the total charge of the sphere, and Φ(d) = C elstat d −1 is the well-known Coulomb potential with its inverse distance dependency and constant prefactor C elstat . Following the theoretical considerations in [39], the section-sphere interaction potential lawπ section-sphere is expressed solely by the scalar separation of the section and sphere midpoint positions. According to the detailed study of the accuracy of this simple resulting interaction law in [39], the intentional neglect of the orientation dependence ψ b-s turned out to be a reasonable approximation in the case of circular, homogeneous cross-sections and long-range interactions as considered here. The required variation of Equation (2.3) for obtaining the corresponding virtual work contribution, the subsequent spatial discretization of the beam centerline for arriving at the discrete element force vector for both the beam and the sphere as well as the consistent linearization of these terms will be omitted here for the sake of brevity. They follow directly from substituting the respective expressions in the equations for the SSIP approach as presented in [39].
Note that arguably the most critical limitation of this computational model is the use of a Coulomb interaction potential law, which neglects the presence of counterions in the electrolyte solution and the associated screening of charges that in turn significantly reduces the range of electrostatic interactions in biological hydrogels. This can be seen as a pragmatic, simplified model chosen due to the current lack of a sphere-beam or beam-beam interaction formulation for screened electrostatic interactions, and the implications will be discussed in the following paragraph.

Discussion of neglecting screening effects.
Regarding the impact of this model assumption on the results of this computational study, it is expected that the effect of charges in general and the range of interaction in particular is overestimated and that applying the simple Coulomb interaction model is thus inadequate for predicting the particle mobility in a quantitatively correct manner. However, the qualitative behavior of the particle and the trends in the statistical quantities of interest such as the mean squared displacement (MSD) of the particle over time for varying charge density is expected to be meaningful and thus allow for both valuable insights in the biophysical system behavior and mechanisms as well as a first proof of principle for this novel computational modeling approach in general. Based on the experimental observations and in anticipation of the obtained simulation results, this reasoning is supported by the fact that the most effective trapping mechanism is the one of a persistent, strongly adhesive contact between the particle and the oppositely charged network fibers, which is a scenario with minimal surface separation and thus minimal screening effect. Therefore, the behavior of a particle with medium to large distance to the nearest fibers is thus expected not to be reproduced correctly by the Coulomb interaction model, whereas in contrast, the practically much more important regime of small separations should be represented with sufficient accuracy in order to allow for the aforementioned analysis of trends and basic mechanisms. In order to still account for screening charges, the cutoff distance of the interaction is set to r cutoff = 2 µm, which is defined via the separation of the sphere and fiber midpoint position and thus effectively neglects any interaction forces beyond particle-fiber surface separations of g cutoff ≈ 0.96 µm. Finally, despite the fact that the development of a screened electrostatic interaction formulation for instance based on the Debye-Hückel approximation of the Poisson-Boltzmann theory would go beyond the scope of this computational study, it is clearly considered an important and promising future extension of the present model that should be used for both subsequent verification of the drawn conclusions and specific analysis of the influence of salt concentration as well as ion-specific effects observed in experiments [12].
The parametrization of the electrostatic interaction model used throughout this study is given as follows. Based on assuming the dielectric permittivity of water (at room temperature) for the surrounding fluid, the constant prefactor of the Coulomb interaction potential law is obtained as C elstat ≈ 1.12 × 10 2 aJ µm fC −2 . As a first step and in accordance with all previous computational studies, the surface charge density of the fibers as well as that of the particle is assumed to be homogeneous and constant along the fibers. In view of the complex, inhomogeneous molecular architecture and thus charge distribution of individual biopolymer filaments and moreover the complex constitution of a real biological hydrogel, this model assumption is once again expected to have a significant influence on the quantitative accuracy of the results. However, it should still allow for qualitative analysis of trends and mechanisms, as argued above in the context of the electrostatic interaction model. A potential improvement on this point is rather a question of detailed experimental fiber characterization and model parametrization than method development, because the sphere-beam interaction model proposed above is capable of describing varying line charge distributions along the filaments. As a first step, however, a constant, homogeneous line charge density of the fibers λ f = −0.25 fC µm −1 is assumed throughout this article and the positive surface charge of the particle will be varied in Section (b) to study its influence on the particle mobility. Finally, a total of 10 integration points per beam element are used to evaluate the contributions of the electrostatic line force distribution along the fibers by means of Gaussian quadrature.

(ii) Repulsive steric interactions
The line-to-line contact formulation proposed in [54] effectively precludes any noticeable penetration of fibers for arbitrary mutual orientations and deformations. It thus serves as the starting point for the dimensionally reduced case of beam-to-rigid sphere contact, which is outlined as follows: Postulating a beam-sphere penalty force law as a linear function of the minimal surface-surface separation (i.e. gap) g b-s with constant scalar prefactor (i.e. penalty parameter) ε b-s yields the two-body interaction potential: Here, the crucial difference to the beam-beam scenario lies in the way the gap g is computed.
Whereas a unilateral closest-point projection is required in the beam-beam scenario, the problem simplifies significantly in the case of a rigid sphere, because the minimal surface separation between the current beam cross-section and the sphere may be expressed in good approximation as: As in the previous section, the required variation of Equation (2.5) for obtaining the corresponding virtual work contribution and subsequent spatial discretization of the beam centerline for arriving at the discrete element force vector for both the beam and the sphere as well as the consistent linearization of these terms follows in close analogy to the beam-beam scenario and will therefore not be presented here.
The parametrization of this contact model is chosen as follows: Throughout this study, a constant penalty parameter ε b-s = 100 pNµm −1 is chosen, which turned out to be sufficiently large such that the maximum penetration of particle and fiber is smaller than 5% of the fiber diameter D f even in the most challenging scenarios of strongly adhesive contact and sudden stochastic forces on the particle in the direction toward the fiber. Finally, 15 collocation points per beam element are used for the evaluation of the contact line force distribution along the fibers.

Remark on fiber-fiber interactions.
At the end of this section, note that contact and electrostatic interactions between fibers are not considered here because it turns out that their mutual separations and orientations are almost entirely determined by the rigid connections of their endpoints at the network vertices. As mentioned earlier, beam-beam interaction formulations for both steric repulsive and electrostatic interactions are readily available and could be directly added as a future extension of this model, however, at the cost of an increased computational effort.

(e) Temporal discretization
In addition to the already discussed spatial discretization of the fibers via beam finite elements, the problem is discretized in time via an implicit direct time stepping scheme. Specifically, the resulting system of first-order stochastic 3 partial differential equations in time is discretized by means of a Backward Euler scheme, as proposed in [30]. Starting from a default time step size of ∆t = 10 −3 s, an adaptive time stepping scheme is applied, which is especially important for resolving the highly nonlinear dynamics during (adhesive) contact interactions and potentially large sudden changes in the magnitude and direction of the thermal forces. The total simulation time per run is set to t end = 20 s, which generally leads to a required number of time steps in the range of n step ∈ [2 × 10 5 , 3 × 10 5 ]. In order to account for the stochastic nature of the thermal forces driving the particle motion, generally two or more random realizations for each -otherwise identical -set of input parameters (i.e., identical, random network geometry and interaction parameters) are computed and considered in the analyses presented in Section 3. 3 Recall that our Brownian dynamics model includes the stochastic thermal forces acting on the spherical particle as stated in

(f) Boundary conditions
As already mentioned in the context of the network modeling (Section (b)), the concept of a representative volume element is used to limit the influence of artificial boundary effects. For this purpose, periodic boundary conditions are applied at each side of the simulation box shown in Figure 3(a). Thus, once the particle as well as (parts of) the fibers leave the representative volume element at any side, they reenter it on the opposite side. Moreover, the steric as well as electrostatic interactions are also evaluated across periodic boundaries. In the majority of the simulations considered in Section 3, no other boundary conditions are applied. However, in a certain batch of simulations, the entire network of fibers will be fixed by means of Dirichlet boundary conditions in order to serve as a reference solution mimicking the limiting case of rigid fibers.

(g) Postprocessing and quantities of interest
Since this study focuses on the (hindered) diffusive mobility of particles, the most important raw data obtained from the simulations is the particle midpoint position in every time step. From this point on, the postprocessing procedure is equivalent to experiments that track the motion of individual fluorescent tracer particles (e.g., [7,8,12]). Based on the discrete time sequence of particle positions, the mean squared displacement (MSD) < ∆r 2 p (τ ) > is computed for any desired time interval (that can be observed in a given simulation run) τ ∈ [∆t, t end ] as follows: Here, Nτ denotes the number of all distinct (but possibly overlapping) time intervals τ obtained for one simulation run. Given that the number of independent samples obtained for large time intervals is naturally limited, only the first 10% of the maximal possible time intervals, i.e., τ ∈ [∆t, 0.1 · t end ], will be considered in the statistical analyses. However, the remaining data is included and indicated by a gray background in all the MSD plots to be presented. Moreover, the mean and standard deviation of the MSD obtained for several realizations of the random network geometry as well as several realizations of the stochastic process of thermal forces will be considered. As the majority of the considered scenarios shows a subdiffusive behavior, the MSD curves over the time interval will be presented and discussed instead of the (apparent) diffusion constant, which obviously depends on the considered time interval and could still be computed from the MSD curves if desired. Other simulation results, such as the resulting axial strains of the fibers, will be presented and discussed for a few specific investigations wherever needed for interpreting the system behavior.

Results and discussion
The following discussion of simulation results is divided into three parts. First, in Section (a), we study the influence of collisions between the particle and the fiber network on the particle mobility. Second, the effect of additional attractive electrostatic interactions between the particle and the fiber network will be investigated in Section (b). Lastly, Section (c) analyzes the special role of fiber stiffness in the case of the most effective hindrance mechanism observed in the simulations.

(a) The effect of solely repulsive steric interactions
To begin with, the double-logarithmic plot of the particle's MSD < ∆r 2 p > as a function of the time interval τ shown in Figure 4   standard deviation of the mean over ten random realizations as indicated by the error bars is negligible within the first 10% of the time interval range (indicated by the white background) that is considered in the following analyses.
Steric hindrance is insignificant as long as particle diameters are smaller than the smallest mesh sizes.
Turning to the influence of repulsive steric interactions, i.e., collisions between the particle and the fiber network, Figure 4(a) shows that a medium fiber volume fractionV f ≈ 0.4% (cyan lines, see Figure 1(a) for an example of the network) has almost no perceptible influence on the MSD.
Only above approximately τ = 0.5 s is a very subtle subdiffusive behavior observable for the three individual realizations plotted here, which indicates the occurrence of a few collisions if the particle travels over longer periods of time. On the one hand, this behavior of almost free diffusion of the particle is expected from the range of cell diameters 2. irregular polygonal network compared to the particle diameter Dp = 1 µm. Moreover, previous experimental studies have made very similar observations for almost neutral particles or high salt concentrations that effectively shield any electrostatic interactions [7,8,12]. However, this is the first computational study with a realistic, irregular polygonal fiber network geometry and thus an important confirmation that the effect of solely steric hindrance is indeed negligible in this regime, where the ratio of network mesh size(s) and particle diameter is greater than one.
Particles with diameters in the range of the mesh sizes are caged in polygonal compartments of random size and show confined diffusive behavior. To push this to the limit where purely steric hindrance and thus filtering by particle size will become significant, a ten times larger number of Voronoi points n VP = 600 corresponding to a fiber volume fraction ofV f ≈ 1.8% and cell diameters of 1.0 − 2.5 µm, i.e., in the order of the particle diameter Dp = 1 µm, have been applied. An example of the resulting network architecture is illustrated in Figure 3(a). To probe the most extreme effect of steric hindrance, the fibers are completely fixed in space for these two sets of three random realizations each.
The results are plotted in Figure 4(a) (green dotted lines). One of the sets (I) shows close to normal diffusive behavior on very small time scales up to approximately τ = 0.05 s and eventually reaches a plateau value of < ∆r 2 p >≈ 0.7 µm 2 beyond τ ≈ 1 s. The other set (II), which are random realizations using the identical fiber network geometry, but a different initial position of the particle, shows a significantly subdiffusive behavior already for the smallest considered time intervals and quickly reaches a plateau value of < ∆r 2 p >≈ 0.05 µm 2 for any time interval longer than τ ≈ 0.1 s. This is an expected behavior for the diffusion of particles in an irregular network with cell diameters of the same order as the particle diameter and therefore randomly connected sufficiently large cells that together form a polygonal volume surrounding the initial position of the particle, which the particle cannot leave under any circumstances. The specific compartments that the particle is able to explore starting from either of the two initial positions in the identical network are illustrated in Figure 4(b) and 4(c). Here, the overlay of all particle positions throughout the entire simulation is shown in dark blue for the case of the higher MSD plateau value (I) and in orange for the case of the lower MSD plateau value (II) observed in Figure 4(a). This behavior is known as confined diffusion and has been theoretically described and experimentally observed e.g., in the context of studying cadherin molecule mobility in plasma membranes [55]. As discussed for instance in [1], such a filtering mechanism based on size clearly has a very effective selectivity and is applied by organisms to strictly preclude the access of any objects larger than the characteristic mesh size, e.g., of the nasal mucous membrane.
However, there are still open questions concerning the mobility of (medium to) large objects in biological hydrogels taking into account the continuous reorganization of biopolymer networks based on both (de-)polymerization and the transient nature of crosslinks. Such a transport mechanism for relatively large particles has been observed in experiments [6] and recently been investigated also in a theoretical and computational model [35], where crosslink binding dynamics are influenced by the diffusing particle. It is also suggested that this kind of mechanism could play a role in the selective permeability of the nuclear pore complex, for which the governing principles are still under debate (see e.g., the review article [1]). Replacing the Voronoi-type network in the present computational model by that of a transiently cross-linked, self-assembled network driven by Brownian motion [32,33] is thus considered a promising future step.

(b) The effect of additional attractive electrostatic interactions
In addition to repulsive steric interactions, attractive electrostatic interactions due to uniformly distributed, opposite charges on the particle and the fibers will now be considered. This has been confirmed to be the most effective hindrance mechanism for particles smaller than the mesh sizes both in experiments (e.g., [7,8,12]) and simulations (e.g., [27,28]). Throughout this section, we thus keep the fiber volume fraction fixed at the medium valueV f ≈ 0.4%, which has been shown to not have a noticeable influence on the particle's MSD in Figure 4(a) (cyan lines).

On average, the degree of subdiffusion increases with the strength of attraction.
The resulting MSD curves for a low (green), medium (cyan), and high (red) value of the particle's surface charge Qp are shown in Figure 5 and compared to the analytical reference solution for free diffusion (black dashed line). As mentioned earlier, the fiber volume fractionV f ≈ 0.4% is identical for all simulation runs. In particular, 5 different random network geometries with 2 random realizations of the stochastic forces each have been simulated for each of the three different charge values. Figure 5(b) shows these 10 independent realizations for each particle charge value and the corresponding mean values and standard deviations are plotted in Figure 5(a). On average, the degree of subdiffusion increases with the strength of attractive interaction, which has been suggested by previous experiments (e.g., [7,12]), and has been confirmed by previous computational studies using ordered (e.g., [26,27]) and unordered [28] arrays of straight, rigid, mutually orthogonal fibers. Here, the smallest charge value Qp = 0.125 × 10 −3 fC leads to a very small degree of subdiffusion, which is in fact quite similar to the one observed in the limit of no charge shown in Figure 4(a). In contrast, medium (cyan) and high charge values (red), which are a factor of 8 and 64 higher than the smallest charge value, significantly hinder the diffusion of the particles already on (very) small time scales τ < 0.1 s. In this regime, the slopes of the MSD curves are however close to one, which suggests normal diffusive behavior with a decreased diffusion constant. slopes in the broad range from zero to one, and even some examples for superdiffusive behavior with a slope greater than one can be observed for medium charge values. Such a significant increase in the variation of MSD values as well as slopes is in excellent agreement with the experimental results from [12]. Their subsequent analysis of squared displacement values over time revealed the existence of sudden trapping and escape events of a particle switching between an almost free, a loosely bound, and a tightly bound state. This hypothesis is confirmed by the simulation results shown in Figure 6(a), where very similar motion patterns can be observed. (Strongly) charged particles jump between local aggregates of fibers, i.e., opposite charges. While a causal link of trapped states and local aggregation of fibers, i.e., charge patches of opposite sign, has already been suggested in several experimental studies (e.g., [7,12]), the limited spatio-temporal resolution of single particle tracking and imaging has so far precluded a direct proof. The particle trajectory shown in Figure 6(b) and the detail view of the fiber aggregates corresponding to the three observed trapped states illustrated in Figure 6(c) clearly confirm this causality. Related observations have been made in the recently published computational study [28], which for the first time considered disorder in the still rigid, straight, mutually orthogonal fibers and described a resulting, so-called dense-region trapping. In the present study, the random, irregular polygonal network structure, very similar to the one of real biological hydrogels (cf. Figures 1 and 3), allows for a more detailed analysis of this trapping mechanism. Particularly, these fiber aggregates are found at network vertices and small polygonal faces in the network that have a lateral spatial extension smaller than the particle diameter. Figure 6 also suggests that the magnitude of the remaining thermal fluctuations observed in the squared displacement of the particle is a measure for the strength of the trapping and that this strength is proportional to the local fiber, i.e., charge density. Finally, the results suggest that the stronger the particle is immobilized, i.e., the smaller the fluctuations in the squared displacement are, the longer the particle will remain at this location.

The variability of MSD values and slopes increases for longer time intervals, which indicates that particles randomly switch between distinct motion patterns.
The irregularity in the fiber/charge distribution gives rise to three distinct motion patterns, one of which allows the particle to travel via successive jumps. From a mechanical point of view, the spatially varying distribution of fibers and thus charges gives rise to a random 3D potential field to be explored by the particle (similar to the one used in our previous contribution [31] to study the effect of filament prestress in biopolymer networks). The location and values of the (local) minima in the potential field as well as the potential barrier and paths connecting them strongly depend on the characteristic geometrical properties of the network such as connectivity, distribution of fiber segment lengths and the resulting distribution of cell and mesh sizes. Considering real biological hydrogels further extends the list of crucial influencing factors to the type and fractions of load-carrying components and their specific surface charge distributions. All of these factors that shape the characteristic potential landscape will determine whether and how the stochastic thermal excitation causes the particle to either (1) remain at one location being completely immobilized, (2) cycle between neighboring minima being restricted to a certain region, or (3) travel through the hydrogel via jumps -potentially also over long distances. All three motion patterns (1)-(3) can be identified among the realizations shown in Figure 5(b). For long time intervals τ > 0.1 s, both patterns (1) and (2) lead to a similar behavior of confined diffusion as obtained for the case of caged particles considered in the previous section 4 . The plateau MSD value thus allows to draw conclusions with respect to the volume enclosed by the particle's initial position and the location(s) of the potential minimum/minima (i.e., fiber/charge agglomerations) that the particle has visited. By looking at the sequence of snapshots over the entire simulation time, it has been verified that those four (highly charged) particles with the smallest MSD plateau values of < ∆r 2 p >≈ 2 × 10 −2 µm 2 are completely immobilized at one single location in the network (cf. motion pattern (1)). In contrast, those medium and highly charged particles with intermediate values for the MSD slopes jump or smoothly transit between local fiber agglomerations in cycles (cf. motion pattern (2)) such as observed in the example shown in Figure 7. Finally, also motion pattern (3) has been identified in the set of individual realizations, as already shown in Figure 6 for a particle with high surface charge. This case can easily be identified among the MSD curves plotted in Figure 5(b) as the one with the largest MSD values for long time intervals.
To conclude this section it can be stated that strong attractive forces mostly lead to motion patterns (1) and (2) that effectively immobilize the particles in a confined volume. However, the motion pattern (3) allows -at least theoretically -for a travel of particles over considerable distances by means of a series of successive jumps. The effectiveness of this transport mechanism depends on the irregularity of the fiber/charge distribution in space and in particular the (relative) height of the potential barriers between the fiber/charge agglomerations. While the particle in the specific example shown in Figure 6 seems to be completely immobilized after two consecutive jumps, one might think of another special design of the fiber/charge distribution with a certain degree of periodicity that could lead to an effective transport of the particle also  over large distances and in small time intervals. Such a directed motion would maximize the mobility of a diffusing particle.
Altogether, the results of this section indicate that the irregularity of the potential landscape exerts a crucial influence on the overall effectiveness of particle immobilization and therefore on the selectivity of hydrogels. The aforementioned factors that shape the effective potential field explored by the particle are known to vary substantially between the multitude of different biological hydrogels (see e.g., [1]). Therefore, a systematic parametrization (and if required an extension) of the present computational model with respect to other classes of gels is expected to be a valuable means for studying the species-specific variations of the general behavior and principles observed so far. In the long run, this might even lead to simulation-based prediction tools enabling a case-specific choice or design of drug delivery vehicles. Being able to optimize the properties of the carrier to effectively attach the drug to the carrier and yet effectively diffuse through the body would clearly be invaluable in this context.
Brief discussion of computational aspects. In total, 60 simulations with at least 2 × 10 4 time steps each have been conducted and evaluated for the results shown in this section (cf. Figure 5(b) for all realizations). One simulation typically took 2-4 days if run in parallel on 16 cores 5 on a Linux cluster. The main drivers for the computational cost are the fiber volume fraction and the presence and strength of electrostatic interactions leading to a complex interplay of repulsive steric and attractive electrostatic forces that require a fine temporal resolution (i.e., smaller time steps) and make the nonlinear problem more challenging to solve.

(c) The influence of fiber stiffness/compliance
Up to this point, we haven't discussed the influence of the fiber stiffness, which mainly influences the amount of fiber deformation and -besides the realistic network geometry -is the second unique feature of the present computational model. In the problem setup considered here, fiber deformations originate exclusively from contact and/or electrostatic interactions with the Brownian particle. Thus, the highest particle charge Qp = 8 × 10 −3 fC will be considered in this section, because the most frequent and strongest interactions can be expected in this case. To determine the point where fiber deformations begin to change the results, the value for Young's modulus has been varied systematically starting from the theoretical limit of rigid fibers as outlined already in Section (b). In this section, the results for rigid fibers will thus be compared to those obtained for a value of E = 0.1 MPa where the first differences can be observed, and to the ones obtained for a ten times larger value E * = 10 · E = 1 MPa, which are basically identical to the case of rigid fibers. Note that -just as in the previous section -the fiber volume fraction will be kept constant at the medium valueV f ≈ 0.4% throughout this section. Figure 8 compares the MSD curves obtained for the low fiber stiffness resulting from Young's modulus E (red lines) with the ones obtained for a ten times higher value for Young's modulus E * (yellow lines) and the ones obtained for the limit of rigid fibers (cyan dashed lines). Apart from the fiber stiffness, the compared sets of ten realizations each are identical, in particular with respect to the 5 different random network geometries and the 2 different sequences of the random stochastic forces each.

The influence of fiber stiffness on particle mobility is insignificant within the range of stiffness values reported for ECM gels.
There is no perceptible difference between the results for the medium fiber stiffness using E * and those for the limiting case of rigid fibers. Recall from Section (b) that the value E * = 1 MPa is already the lower limit of the wide range E = 1 MPa − 1 GPa of values reported for experiments with collagen gels [48,49].
For more flexible fibers, we observe and speculate about a few competing effects that seem to cause an overall decrease of particle mobility. First, let us look at the resulting overall effect in terms of the change in ensemble-averaged MSD curves. For the lowest fiber stiffness E = 0.1 MPa considered in our simulations, the curves are shifted towards smaller MSD values (see Figure 8). The difference is not significant due to the considerable variability already observed in the results of the last section. However, these results allow to conclude that -for the given set of parameters -a value below E * will begin to influence the simulation results. As an explanation for this behavior, we speculate that the softer fibers . . lead to an increased adhesive contact area because their shape adapts to the surface shape of the particle and that this behavior in turn traps the oppositely charged particles more effectively. In addition, the softer fibers will presumably follow the particle's thermal excitation more closely and thus lead to smaller peak values of particle accelerations and interaction forces that act to separate the adhesive contact between fiber(s) and the particle, such that escape events will become less likely. Both effects might contribute to the observed overall reduction of particle mobility for fibers with a low stiffness. This reasoning is supported by the previously recognized high importance of trapping and escape events as concluded from the results of the previous Section (b).
Generally, there might be more mechanisms how very soft fibers can influence the particle mobility, including some that increase particle mobility. In an effort to shed some light on this, we look at the direct, pair-wise comparison of a pair of simulations with different fiber stiffness values E and E * and an otherwise identical setup, i.e., in particular, with identical network geometry, identical initial particle position, and identical sequence of stochastic forces (see Figure 9). The first and second pair of realizations with otherwise identical parameters (first pair marked with crosses and second pair marked with diamonds) show nearly identical MSD curves for E (red) and E * (yellow), respectively. Conversely, the third pair of realizations with otherwise identical parameters (marked with circles) leads to entirely different results for E (red) and E * (yellow). This leads to the conclusion that small differences in the fiber behavior can trigger entirely different particle trajectories in individual realizations of the stochastic process. Given the already large variability in the results of individual random realizations for identical fiber stiffness observed in the last section (see the right half of Figure 8(b)), this sensitivity is not surprising. The stochastic nature of the problem impedes the direct investigation of differences in particle behavior by means of a pair-wise comparison and underlines the importance of the already analyzed ensemble-averaged results shown in Figure 8(a) and discussed above. However, by taking a closer look at the second pair of realizations and watching the two resulting particle motions as overlaying videos, we can still investigate a quite obvious mechanism that influences the particle mobility providing, in this case, nearly identical MSD curves. Indeed, one would intuitively argue that stiffer fibers constrain the motion of a particle more effectively as long as it remains in the state of adhesive contact with the fiber, because the softer fibers deform due to the thermal forces acting on the trapped particle. This situation is shown in Figures 10(b) and 10(a), which show an overlay of fiber configurations over all time steps for these two realizations. The two corresponding MSD curves in Figure 9 indeed confirm that the softer fibers (red diamonds) allow for slightly higher MSD values than the stiffer fibers (yellow diamonds) on time intervals τ ≥ 0.1 s. However, these MSD curves also show that this effect is rather negligible in terms of overall mobility of the particles. This is also supported by recognizing the small magnitude and localized extent of axial strains in the network (at an exemplarily chosen point in time) in Figure 10(c).
To conclude this first brief study of how deformations of semiflexible fibers influence the mobility of oppositely charged particles, it can thus be stated that our results indicate a -rather counter-intuitive -overall decrease in particle mobility if compared to (almost) rigid fibers. We suggest that this is the result of a decreased probability for escape events due to a) an increase in the adhesive contact area between particle and fiber and b) a higher ability of the fibers to follow the thermal excitation of a trapped particle, thus reducing the peak values of particle accelerations and interaction forces that act to break the adhesive binding. It is important to note, however, that this influence of the fiber stiffness on the particle mobility is negligible in the range of stiffness values that have been reported for the collagen I fibers prevailing in ECM gels. Within this range, the fibers do not show any noticeable deformations for the scenario considered here and the associated influence on the particle trajectory thus vanishes. Nevertheless, the influence of fiber deformations is likely to be observed for other types of (biological) hydrogels with very thin or soft and therefore more flexible filaments such as mucin or F-actin, and maybe even for ECM gels as a result of dysregulated fiber stiffness. A more detailed investigation of this aspect, both in simulations and experiments, is thus considered a promising avenue of future research.

Conclusions and outlook
This article presents the first computational study of the diffusive mobility of particles in hydrogels with a realistic fiber network model. It proposes a novel computational approach based on, most notably, the modeling of individual, deformable fibers via the beam theory, the Voronoi tessellation of the periodic simulation box to obtain random, irregular network geometries, and the beam-sphere interaction model for contact and electrostatic interactions. In the case of only repulsive steric interactions, it is found that the hindrance of the particle mobility is insignificant as long as the mesh sizes of the fiber network are larger than the particle diameter. If, however, the mesh sizes are in the order of the particle diameter, the particle is effectively caged in a polygonal fiber hull of random shape and size and shows a behavior known as confined diffusion and is characterized by a plateau in the mean squared displacement (MSD) curve for long time intervals. Within the given problem setup and the focus on relatively small particles, these are expected results validating the novel approach. However, this effect of steric hindrance will become highly relevant if the effective transport of relatively large particles as observed in experiments [6] is considered. In this context, including the dynamics of the fiber network such as its self-assembly driven by the Brownian motion and transient reorganization of cross-links (as demonstrated for the directly compatible computational model applied in [32,33]) might be an important model component as suggested by recent findings in the context of the nuclear pore complex [9,35] or the dynamic secretion and shedding of mucus layers [36].
Turning to the effect of additional electrostatic interactions between the fibers and the oppositely charged particle, the prevailing notion that the degree of hindrance on average increases with the strength of attraction has been confirmed by the numerical experiments with five different network geometries and two random realizations each. Moreover, an increased variability of the particle's mean squared displacement values and slopes in the regime of long time intervals has been observed and excellently agrees with previous experimental results [12]. A detailed look at the 3D particle trajectories within the fiber network provides a first direct proof for the existence of distinct motion patterns of the particles, which explains the variability in the MSD curves. As hypothesized in the previous work [12], the particles stick to oppositely charged fiber/charge aggregations experiencing more or less strong trapping and eventually escape due to the ongoing thermal excitation, only to be quickly attracted to another fiber/charge aggregation. While some particles remain completely immobilized at one and the same location for the entire 20s of simulation time, others smoothly or rapidly cycle between two local minima in the potential landscape. Both of these motion patterns lead to a behavior on longer time intervals that is very similar to the confined diffusion for caged, uncharged particles as described above. However, the diffusive mobility on short time intervals is significantly reduced as well due to sticking to the fibers. The third motion pattern observed is the one of several successive jumps that -at least theoretically -could serve as a transport mechanism also over longer distances if the potential landscape is formed accordingly and e.g., shows some degree of periodicity and directional preference.
Altogether, these findings indicate that the precise shape of the effective 3D potential field explored by the particle has a crucial influence on its mobility. In view of the broad variety of biopolymer hydrogels with diverse chemical compositions and biophysical properties, the current computational model could thus be leveraged to study the individual selective filtering behavior for a large number of particle-hydrogel property combinations. Based on the recognized importance of the precise fiber/charge distribution in the system, two points seem to be of particular importance to achieve a case-specific, highly accurate and reliable prediction. First, the inhomogeneous charge distribution along the fiber should be both determined (e.g., by experimentally analyzing the molecular architecture) and applied in the model. Second, the specific composition and geometry of fiber networks should be determined (e.g., by processing electron microscopy images) and applied in the model. Also the inclusion of the dynamic selfassembly and reorganization of networks mentioned above would be worth considering in this respect.
As a last particular aspect investigated in this study, fiber deformations have been found to be negligible within the range of realistic values for the stiffness of collagen I fibers prevailing in ECM gels. To be more precise, varying the value for the Young's modulus over the broad range of reported values for the considered ECM gels has led to identical results as obtained for the theoretical limit of rigid fibers. If, however, more flexible fibers are considered, our simulation results indicate an overall decrease of particle mobility if compared to (almost) rigid fibers - an outcome that is rather counter-intuitive. We suggest that this is the result of a decreased probability for escape events due to a) an increase in the adhesive contact area between particle and fiber and b) a higher ability of the fibers to follow the thermal excitation of a trapped particle, thus reducing the peak values of particle accelerations and interaction forces that act to break the adhesive binding. In real systems, these trends might be observed for other types of (biological) hydrogels with very thin or soft and therefore more flexible filaments such as mucin or F-actin, and maybe even for ECM gels as a result of dysregulated fiber stiffness. A more detailed investigation of this aspect, both in simulations and experiments, is thus considered a promising avenue of future research.
In addition to the presented simulation results and the gained insights, this study provides an extensive proof of concept for the application of the novel computational model. As outlined above, in the short to medium term many important findings especially for various particle/hydrogel-specific behaviors and mechanisms can be expected from applying different parametrizations and, additionally, from integrating the suggested model extensions. In the long term, further validation and advances of the present computational model toward a truly predictive tool could ultimately lead to a case-and patient-specific choice or even design of pharmaceuticals and also to a case-and patient-specific assessment of infection risk.
Ethics. This article does not contain any studies with human or animal subjects.
Data Accessibility. All required information to reproduce the results of the numerical experiments is provided or referenced in this article. The source code is part of a wider code base [41] that will be made available in the future. The data is provided as electronic supplementary material.