A Novel Modeling and Simulation Approach for the Hindered Mobility of Charged Particles in Biological Hydrogels
Maximilian J. Grill, Jonas F. Eichinger, Jonas Koban, Christoph Meier, Oliver Lieleg, Wolfgang A. Wall
rrspa.royalsocietypublishing.org
Research
Article submitted to journal
Subject Areas: computational biophysics, modelingand simulation, biomaterials
Keywords: hindered particle diffusion, biologicalhydrogels, electrostatic interaction,deformable fiber network, beamtheory, finite element method
Author for correspondence:
Maximilian J. Grille-mail: [email protected]
A Novel Modeling andSimulation Approach for theHindered Mobility of ChargedParticles in BiologicalHydrogels
Maximilian J. Grill , Jonas F. Eichinger ,Jonas Koban , Christoph Meier , OliverLieleg and Wolfgang A. Wall Institute for Computational Mechanics, TechnicalUniversity of Munich, Germany Munich School of Bioengineering, TechnicalUniversity of Munich, Germany
This article presents a novel computational model tostudy the selective filtering of biological hydrogelsdue to the surface charge and size of diffusingparticles. It is the first model that includes therandom 3D fiber orientation and connectivity of thebiopolymer network and that accounts for elasticdeformations of the fibers by means of beam theory.As a key component of the model, novel formulationsare proposed both for the electrostatic and repulsivesteric interactions between a spherical particle and abeam. In addition to providing a thorough validationof the model, the presented computational studiesyield new insights into the underlying mechanismsof hindered particle mobility, especially regardingthe influence of the aforementioned aspects that areunique to this model. It is found that the precisedistribution of fiber and thus charge agglomerationsin the network have a crucial influence on themobility of oppositely charged particles and givesrise to distinct motion patterns. Considering thehigh practical significance for instance with respectto targeted drug release or infection defense, theprovided proof of concept motivates further advancesof the model toward a truly predictive computationaltool that allows a case- and patient-specific assessmentfor real (biological) systems. © The Author(s) Published by the Royal Society. All rights reserved. a r X i v : . [ phy s i c s . b i o - ph ] J a n r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
1. Introduction
The remarkable ability of hydrogel forming biopolymer networks to control the mobility ofdifferent kinds of diffusing molecules and particles individually is of crucial importance fornumerous functions of the human body. See [1] for a recent review of this topic. On the onehand, this selective permeability gives rise to the protection of the organism against pathogenssuch as viruses that are effectively hindered from invading and traversing the organism. On theother hand, it ensures the effective transport of a broad variety of substances that are usefuland important for the organism. Examples of such biological hydrogels include mucus, theextracellular matrix (ECM), intracellular biopolymer networks (comprising, e.g., neurofilamentsor actin), the vitreous humor and the matrix of biofilms, and can thus be found throughout theentire human body and in countless other places in nature. The high practical relevance of thisfundamental topic extends over multiple fields ranging from medical diagnosis to the therapyof body malfunctions and targeted drug release. This creates yet unimagined possibilities intechnical 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 bysize – the surface properties of the particles also plays an important role (see e.g., [1–16]). Thisnew paradigm thus suggests filtering not (only) by size but also through a combination of otherparticle properties such as charge, hydrophobicity or binding affinity and can thus be referred toas interaction filtering. Despite considerable scientific effort in this field, many aspects concerningthe underlying mechanisms and specific conditions remain unknown. To a large extent, thiscan be explained by the following three factors. First, there is a great complexity in the manydifferent biological systems both in number and diversity of components, e.g., with respect totheir molecular architecture. Second, there are big challenges and limitations with respect toexperimental preparation and measurement techniques when it comes to the required high spatialand 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 onthe molecular scale. This lack of detailed, fundamental microscopic understanding prevents anyreliable prediction of the diffusive mobility of molecules and particles in other than the fewparticle-hydrogel systems already studied in vitro. At this point, the development of accurateand efficient computational models, which are capable of resolving small scales and coveringlarge spans in space and time, is expected to substantially contribute to scientific progress in thisfield.Compared to the large body of literature reporting on experimental work in this area, relativelyfew computational studies have been published so far. Most of the early and also some ofthe recent work focused solely on the excluded volume effect of the fibers, i.e., the repulsivesteric interactions with the network fibers that hindered the free diffusion of a (hard) sphericalparticle (see e.g., [17–22]). In his Monte Carlo simulations, Saxton [23] was the first to includeand study other than steric interactions in the form of a binding model. Since the recognitionof the dominant role of electrostatic and possibly other types of molecular interactions asoutlined above, several computational studies have included electrostatic effects and confirmedthe trends observed in experiments and shed light on the underlying mechanisms [13,24–28].Particularly the recent works published by Hansing et al. [27] and Hansing and Netz [28] werevery 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 foundgood agreement of the simulation results with several sets of experimental data, which confirmsthe validity of such modeling approach.The computational model proposed in this work aims to improve especially the modeling ofthe 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 r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... straight rigid fibers [13,27] has been assumed. Zhou and Chen [24] likewise applied a cubic latticeconsisting of beads placed at the vertices and connected by linear spring elements. Hansing andNetz [28] were the first to break the strong geometrical order in all of these models, they onlyallowed, however, for a variation of the spacing of the still straight and infinitely long rigid fiberswith a mutually orthogonal orientation. Yet they found a significant and fundamental differencein the trapping mechanism of ordered and disordered fiber lattices, which can be attributed tolocally denser regions of the network in the case of attractive particle-network interactions. Thisis a strong motivation to work toward a more realistic modeling of the fiber network as a crucialpart 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 beamtheory, thus allowing for arbitrarily curved initial shapes of the fibers and possibly largedeformations, including all six modes of axial strain, (2x) shear, torsion and (2x) bending. Inaddition, the initial spatial distribution, orientation and interconnection of fibers will be modeledas a random 3D Voronoi network, thus mimicking several important geometrical features of realbiopolymer networks such as a random, spatially variable mesh size distribution, a random fiberlength distribution, random mutual orientations of the fibers, and arbitrary connectivity betweenthe fibers. Several of these attributes are expected to play an important role when it comes to bothpurely steric interactions and the combination with electrostatic interactions and will be studied inSection 3. A similar approach to network generation based on Voronoi tessellation has previouslybeen applied in a number of publications, e.g., to study cell-cell communication in a 2D networkof linearly elastic springs [29].Applying such a sophisticated model to the individual fibers and the network they constitutecomes at the cost of an increased complexity and size of the system of equations to be solved inthe simulations. However, the same modeling strategy based on geometrically exact beam finiteelements describing the biopolymer fibers has previously been applied to model the Browniandynamics of individual semiflexible filaments [30,31] and has been proven to be accurate and alsoefficient enough to study large-scale problems such as the process of self-assembly of (differentmorphologies of) transiently cross-linked biopolymer networks [32,33] as well as their (high- andlow-frequency) rheology [34]. On the one hand, this confirms that the novel modeling approachis suited to study the hindered mobility of particles in hydrogels and on the other hand, thisalready 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 firststep in the present study. Moreover, the dynamics of the network, including the reorganizationof the transient cross-links could be included and would be highly interesting since there isexperimental evidence that particles larger than the mesh size can still diffuse through thenetwork by locally breaking inter-fiber links [6]. This has also been confirmed to be an effectivetransport mechanism in a first computational model [35]. Another example for the dynamics ofhydrogel 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 numericalmethods required to solve it as well as the code framework substantially differs from the onesused by the previous computational studies listed above. This is an inevitable consequenceof the aforementioned modeling of individual biopolymer fibers as elastically deformable,geometrically exact beams. In the present model, a set of well-established numerical formulationsfor beams [37,38] is combined with novel beam-sphere interaction models, which are derivedas 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 betweencurved slender fibers undergoing large deformations and is thus an important prerequisite forthe 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 thecomputational model and the required numerical methods. After the presentation and discussionof the results in Section 3, this article will be concluded with a summary of the findings and anoutlook to promising aspects of future research in Section 4. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
2. Computational model and methods
Given the broad variety of biological hydrogels mentioned above, we chose the ECM gel(s) usedin the comprehensive experimental studies of particle mobility by Lieleg et al. [7] and Arends etal. [12] to serve as the main reference for the specific setup and parametrization of the versatilecomputational model. For these gels, an additional, subsequent characterization in terms oftheir 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 modelconsiders only the following three key components, which are thought to mimic the crucialinfluences of the studied problem: a fiber network, a diffusing particle, and their mutualinteractions. All three model components will be considered individually in the following andselected further aspects of the simulation setup such as boundary conditions and post-processingof the results will be discussed. The software package used for the simulations in this article is theparallel, 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 commonlyused in nonlinear finite element frameworks for structural dynamics. In short, this approachconsists of the following steps: According to the principle of virtual work, the weak form of themechanical balance equations is derived and subsequently discretized in space and time. Given aproper set of initial and boundary conditions, an implicit load/time stepping scheme is appliedand in every step the solution of the resulting discrete system of nonlinear equations is founditeratively by means of Newton’s method. Refer to Ref. [39] for a discussion of selected aspects ofthe 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 3DVoronoi tessellation of the cubic simulation box (size × × µ m ) generated via the opensource 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 ofso-called Voronoi points . The output are the vertices and edges of a random, irregular, polygonalnetwork that are used to define the position and orientation of the initially straight beam segmentsas well as their interconnections. For Voronoi-based tessellation, the connectivity number, i.e., thenumber of fibers branching off at each junction point, is 4, which agrees well with values of 3 to4 reported for ECM gels [44]. In the future, this parameter could also be adapted by randomlyremoving or adding fibers to match other hydrogel architectures. At each junction, the beamendpoint positions and rotations are coupled, which corresponds to a model assumption basedon the expected mechanical behavior of the chemical binding between the fibers, and could beonce again adapted e.g., to hinged connections in a straightforward manner. Also, the bindingand unbinding of such connections at given rates could be included by following the approachdescribed 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 tessellationand thus the resulting network geometry is chosen to be periodic in space. Altogether, a simplevisual comparison with electron microscopy images of real biological hydrogels in Figure 1reveals a high similarity both in terms of the random, irregular, polygonal structure as wellas its characteristic properties, such as the distribution of pore sizes, fiber segment length andconnectivity.The individual fibers are modeled by (geometrically exact) 3D beam theory, assumingundeformable cross-sections of circular shape and a hyperelastic material law. Specifically, we Note that the original term Voronoi particles from [42] is not used here to avoid confusion with the diffusing particle(s). r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... (a)(b)(c) Figure 1.
Network architectures: (a) The result of a random, periodic Voronoi tessellation, which is used as the initial,stress-free configuration in the computational model (box size is × × µ m , rendered using Blender [43]). Theinset shows the scaled system (scale bar indicates µ m ) used to ease the comparison with the other images. (b)–(c)SEM images of four different basal lamina gels used in in vitro experiments (reprinted from [40], scale bars indicate (b) µ m and (c) µ m ). apply the Simo-Reissner beam theory [45–47], which accounts for 6 deformation modes of axialtension, (2x) shear, torsion and (2x) bending. Further specifications such as the dimensional andconstitutive parameters are chosen to mimic collagen I as the key constituent of the targetedclass of ECM hydrogels, however, all these parameters can be easily adapted to study theirinfluences or to model other fiber species. In the present study, the cross-section diameter isset to D f = 75 nm , Poisson’s ratio is set to ν = 0 . , and Young’s modulus is varied from as lowas E = 0 . up to the theoretical limit of rigid fibers to study the influence of the fiber r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... stiffness. Note that experimental measurements for collagen I suggest values in the wide rangeof E = 1 MPa − [48,49], which are covered in this work as well. Based on the resultsin Section (c), it will turn out that the influence of the fiber stiffness on the particle mobilityin the problem setup considered here is negligible for the realistic range of values for Young’smodulus E , and a noticeable difference in results can only be observed below a threshold valueof E ∗ = 1 MPa . It can therefore be assumed that the specific values for the fiber material are ofminor importance and that the deformation of fibers will only become important if the dynamicsof 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 resultingvalues of the mean and standard deviation of the fiber volume fraction ¯ V f and theminimum/average/maximum cell diameters as a function of the number of Voronoi points n VP .For this statistical analysis of the random network geometries, a box size of × × µ m and (a) (b) Figure 2. (a) Mean and standard deviation of the fiber volume fraction (black circles with error bars, which are smaller thanthe symbol size). (b) Mean and standard deviation of the minimum cell diameter (blue pluses with error bars) / averagecell diameter (black circles with error bars, again smaller than the symbol size) / maximum cell diameter (green diamondswith error bars), obtained for three random network geometries for each of the considered numbers of Voronoi points n VP . a fiber diameter D f = 75 nm as stated above is assumed. A number of n VP = 60 Voronoi pointsresult in a fiber volume fraction ¯ V f ≈ . and a range of cell diameters of approx. . − . µ m .This turns out to match the typical mesh size of − µ m reported in [7] quite well and is thuschosen as the default value for most of the simulations conducted in this work (once againrefer to Figure 1 for a visual comparison of model and real hydrogels). The densest networkto be considered in this work is given as n VP = 600 and thus results in a high fiber volumefraction of ¯ V f ≈ . and cell diameters in the range of approx. . − . µ m . This second valueof the network parameter n VP = 600 is motivated by the example of the human amniotic basalmembrane [44], which appears to be much denser than the one considered above. One exampleof the resulting model network is shown in Figure 3.In order to investigate and average out the influence of the specific network geometry, 5different realizations of the random network generation process will be considered as inputfor 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 fibersare assumed to be athermal because the lengths of the individual fiber segments are generallymuch smaller than their persistence length (cid:96) p ≈ µ m (resulting from the parameters for the Here, the cell diameter for each of the n VP cells in the network has been computed as an approximated inscribed sphere ofthe irregular polyhedron based on the shortest distance of the Voronoi point to each of the cell edges. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... (a)(b) (c) Figure 3. (a) Image of an example network resulting for a second, very high value for the fiber volume fraction ¯ V f ≈ . considered in the simulations (box size is × × µ m , rendered using Blender [43]). Comparison of (b) amagnified part of this network with (c) a magnified part of the network with the default fiber volume fraction ¯ V f ≈ . (shown in Figure 1(a)). bending stiffness given above and the thermal energy at room temperature k B T ≈ × − aJ )and thermal undulations will therefore be negligible. However, for another network architecturewith longer fiber segments or a different species of thinner and thus more flexible fibers suchas mucin or F-actin, this effect might become important and may be included in the novelcomputational model by means of the formulation proposed in [30].Finally, the spatial discretization of the fibers makes use of the geometrically exactHermitian Simo-Reissner beam elements proposed in [37], which are mainly based on the well-established element formulation proposed by Crisfield and Jeleni´c [50,51]. The applied centerlineinterpolation using cubic Hermite polynomials ensures both high accuracy in terms of spatialapproximation and a C -continuous geometry representation. This is particularly important forsmooth contact kinematics and smooth interaction force distributions as has been shown in thecontext of beam-beam interactions in [39,52]. When creating the finite element discretization foreach of the random network configurations, a default element length l ele = 1 µ m is used for eachstraight fiber segment and one additional shorter beam element is created for the remainder of therandom segment length if needed. Based on the previous experience (from challenging scenarioswith large deformations such as presented in [39]) with this kind of beam finite element featuring r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... fourth order spatial convergence, this is considered to be a sufficiently fine spatial resolutionfor the magnitude of deformations observed throughout this computational study, such thatthe influence of the spatial approximation error on the results is expected to be negligible. Toconclude, this strategy typically leads to a number of fiber segments n f ≈ , a number of beamelements n ele ≈ , and a number of nodes n node ≈ for n VP = 60 . Likewise, for the densestnetwork with n VP = 600 , we obtain n f ≈ , n ele ≈ , and n node ≈ . (c) Spherical particle The particle is modeled as a rigid sphere and is therefore uniquely described by its midpointposition r p ( t ) ∈ R as a function of the time t as well as its diameter D p = 1 µ m that is chosento comply with the experimental studies presented in [8,12]. As in several of the aforementionedprevious computational studies (e.g., [27,28]), the Brownian dynamics of the particle is modeledby the Langevin equation (see e.g., [53]), including the stochastic thermal as well as viscousdrag forces that are repeated here for the reader’s convenience. The velocity-proportional dragforce f s,visc implicitly models a quiescent surrounding fluid and makes use of the frictioncoefficient of a sphere γ = 3 πηD p according to Stokes: f s,visc = γ I × ˙ r p . (2.1)Here, ˙ r p denotes the particle velocity, and the fluid viscosity is chosen as η = 1 mPa s , whichcorresponds to the viscosity of water at room temperature and has been found to be in goodagreement with experimental tracer studies of freely diffusing particles [12]. Following theconsistent modeling of the Brownian dynamics of slender biopolymers within the space- andtime-discrete theoretical framework of the nonlinear finite element method according to [30], thestochastic thermal forces acting on the sphere are given as: f s,stoch = (cid:112) k B T γ I × ∂ W ( x, t ) ∂ x ∂ t . (2.2)Here, the thermal energy is set to k B T = 4 . × − aJ corresponding to room temperature andthe last term describes the space-time white noise resulting from a two-dimensional Wienerprocess W ( x, t ) . Both contributions from viscous and stochastic forces on the spherical particleare added to the total virtual work of the system, which other than that, includes the contributionsfrom internal elastic forces of the fibers and the contributions from the particle-fiber interactionsthat will be discussed next. Afterward, the temporal discretization and time stepping scheme willbe 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 thespherical 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 modeledat the level of individual particle-fiber pairs and are evaluated as a fully resolved, resultingline force distribution on the fiber. At this point, it becomes clear that the sophisticated fibermodel, including potentially large 3D deformations, which is a unique feature of this novelcomputational model, carries over to interaction modeling in the form of additional challengesto accurately and efficiently describe the sphere-fiber interactions for arbitrarily deformedfiber configurations and mutual orientations. Since the conclusion of previous experimentalas 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 spherevia both repulsive steric and (attractive) electrostatic effects has not been considered in theliterature before. The key idea of our modeling approaches is to consider the sphere-beam r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... interaction as a special, simpler case of the beam-beam interaction. Starting from the formulationfor molecular interactions proposed in [39] and the contact formulation proposed in [54], thefollowing 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 arigid sphere and a beam, are summarized below. Note that hydrophobic interactions are notincluded 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 howeverbe incorporated in a very similar manner as the electrostatic interactions and would be aninteresting 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: Π ia = l (cid:90) ˜ π section-sphere ( r b-s , ψ b-s ) d s, (2.3)where a section-sphere interaction potential law ˜ π section-sphere has been introduced. Such areduced interaction law is an analytical description of the effective interaction of one cross-sectionof 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 spheremidpoint r s 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 interactionpotential ˜ π section-sphere is then numerically integrated along the arbitrarily deformed centerlinecurve of the beam. Here, s ∈ [0 , l ] denotes the arc-length parameter, which is defined in the stress-free, initial configuration of the beam’s centerline curve. This general approach is in close analogyto 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 ofCoulomb interaction. Again in close analogy to the careful choice of the SSIP law for this kind ofinteraction in [39], the following section-sphere interaction potential law ˜ π section-sphere is obtained: ˜ π section-sphere = λ b Q s Φ ( d ) with d = (cid:107) r b − r s (cid:107) . (2.4)Here, λ b denotes the line charge density of the beam, Q s denotes the total charge of the sphere,and Φ ( d ) = C elstat d − is the well-known Coulomb potential with its inverse distance dependencyand constant prefactor C elstat . Following the theoretical considerations in [39], the section-sphereinteraction potential law ˜ π section-sphere is expressed solely by the scalar separation of the sectionand sphere midpoint positions. According to the detailed study of the accuracy of this simpleresulting interaction law in [39], the intentional neglect of the orientation dependence ψ b-s turnedout to be a reasonable approximation in the case of circular, homogeneous cross-sections andlong-range interactions as considered here.The required variation of Equation (2.3) for obtaining the corresponding virtual workcontribution, the subsequent spatial discretization of the beam centerline for arriving at thediscrete element force vector for both the beam and the sphere as well as the consistentlinearization of these terms will be omitted here for the sake of brevity. They follow directly fromsubstituting 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 aCoulomb interaction potential law, which neglects the presence of counterions in the electrolytesolution and the associated screening of charges that in turn significantly reduces the range ofelectrostatic interactions in biological hydrogels. This can be seen as a pragmatic, simplifiedmodel chosen due to the current lack of a sphere-beam or beam-beam interaction formulationfor screened electrostatic interactions, and the implications will be discussed in the followingparagraph. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Discussion of neglecting screening effects.
Regarding the impact of this model assumption on the results of this computational study, itis expected that the effect of charges in general and the range of interaction in particular isoverestimated and that applying the simple Coulomb interaction model is thus inadequate forpredicting the particle mobility in a quantitatively correct manner. However, the qualitativebehavior of the particle and the trends in the statistical quantities of interest such as the meansquared displacement (MSD) of the particle over time for varying charge density is expected tobe meaningful and thus allow for both valuable insights in the biophysical system behavior andmechanisms as well as a first proof of principle for this novel computational modeling approachin general. Based on the experimental observations and in anticipation of the obtained simulationresults, this reasoning is supported by the fact that the most effective trapping mechanism is theone of a persistent, strongly adhesive contact between the particle and the oppositely chargednetwork fibers, which is a scenario with minimal surface separation and thus minimal screeningeffect. Therefore, the behavior of a particle with medium to large distance to the nearest fibersis thus expected not to be reproduced correctly by the Coulomb interaction model, whereas incontrast, the practically much more important regime of small separations should be representedwith sufficient accuracy in order to allow for the aforementioned analysis of trends and basicmechanisms. In order to still account for screening charges, the cutoff distance of the interaction isset to r cutoff = 2 µ m , which is defined via the separation of the sphere and fiber midpoint positionand thus effectively neglects any interaction forces beyond particle-fiber surface separationsof g cutoff ≈ . µ m . Finally, despite the fact that the development of a screened electrostaticinteraction 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 consideredan important and promising future extension of the present model that should be used for bothsubsequent verification of the drawn conclusions and specific analysis of the influence of saltconcentration as well as ion-specific effects observed in experiments [12].The parametrization of the electrostatic interaction model used throughout this study isgiven 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 isobtained as C elstat ≈ . × aJ µ m fC − . As a first step and in accordance with all previouscomputational studies, the surface charge density of the fibers as well as that of the particleis assumed to be homogeneous and constant along the fibers. In view of the complex,inhomogeneous molecular architecture and thus charge distribution of individual biopolymerfilaments and moreover the complex constitution of a real biological hydrogel, this modelassumption is once again expected to have a significant influence on the quantitative accuracyof 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 improvementon this point is rather a question of detailed experimental fiber characterization and modelparametrization than method development, because the sphere-beam interaction model proposedabove is capable of describing varying line charge distributions along the filaments. As a firststep, however, a constant, homogeneous line charge density of the fibers λ f = − .
25 fC µ m − isassumed throughout this article and the positive surface charge of the particle will be varied inSection (b) to study its influence on the particle mobility. Finally, a total of 10 integration pointsper beam element are used to evaluate the contributions of the electrostatic line force distributionalong the fibers by means of Gaussian quadrature. (ii) Repulsive steric interactions The line-to-line contact formulation proposed in [54] effectively precludes any noticeablepenetration of fibers for arbitrary mutual orientations and deformations. It thus serves as thestarting point for the dimensionally reduced case of beam-to-rigid sphere contact, which isoutlined as follows: Postulating a beam-sphere penalty force law as a linear function of the r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... minimal surface-surface separation (i.e. gap) g b-s with constant scalar prefactor (i.e. penaltyparameter) ε b-s yields the two-body interaction potential: Π c ε ,b-s = 12 ε b-s l (cid:90) (cid:104) g b-s ( s ) (cid:105) d s. (2.5)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 problemsimplifies significantly in the case of a rigid sphere, because the minimal surface separationbetween the current beam cross-section and the sphere may be expressed in good approximationas: g b-s ( s ) = (cid:107) r b − r s (cid:107) − R b − R s . (2.6)As in the previous section, the required variation of Equation (2.5) for obtaining thecorresponding virtual work contribution and subsequent spatial discretization of the beamcenterline for arriving at the discrete element force vector for both the beam and the sphere aswell as the consistent linearization of these terms follows in close analogy to the beam-beamscenario and will therefore not be presented here.The parametrization of this contact model is chosen as follows: Throughout this study, aconstant penalty parameter ε b-s = 100 pN µ m − is chosen, which turned out to be sufficientlylarge such that the maximum penetration of particle and fiber is smaller than of the fiberdiameter D f even in the most challenging scenarios of strongly adhesive contact and suddenstochastic forces on the particle in the direction toward the fiber. Finally, 15 collocation points perbeam 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 notconsidered here because it turns out that their mutual separations and orientations are almostentirely determined by the rigid connections of their endpoints at the network vertices. Asmentioned earlier, beam-beam interaction formulations for both steric repulsive and electrostaticinteractions 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, theresulting system of first-order stochastic partial differential equations in time is discretized bymeans of a Backward Euler scheme, as proposed in [30]. Starting from a default time step sizeof ∆t = 10 − s , an adaptive time stepping scheme is applied, which is especially important forresolving the highly nonlinear dynamics during (adhesive) contact interactions and potentiallylarge sudden changes in the magnitude and direction of the thermal forces. The total simulationtime per run is set to t end = 20 s , which generally leads to a required number of time steps inthe range of n step ∈ [2 × , × ] . In order to account for the stochastic nature of the thermalforces driving the particle motion, generally two or more random realizations for each – otherwiseidentical – set of input parameters (i.e., identical, random network geometry and interactionparameters) are computed and considered in the analyses presented in Section 3. Recall that our Brownian dynamics model includes the stochastic thermal forces acting on the spherical particle as stated inEquation 2.2. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... (f) Boundary conditions As already mentioned in the context of the network modeling (Section (b)), the concept of arepresentative volume element is used to limit the influence of artificial boundary effects. Forthis purpose, periodic boundary conditions are applied at each side of the simulation box shownin Figure 3(a). Thus, once the particle as well as (parts of) the fibers leave the representativevolume element at any side, they reenter it on the opposite side. Moreover, the steric as wellas electrostatic interactions are also evaluated across periodic boundaries. In the majority of thesimulations considered in Section 3, no other boundary conditions are applied. However, in acertain batch of simulations, the entire network of fibers will be fixed by means of Dirichletboundary conditions in order to serve as a reference solution mimicking the limiting case of rigidfibers. (g) Postprocessing and quantities of interest
Since this study focuses on the (hindered) diffusive mobility of particles, the most important rawdata obtained from the simulations is the particle midpoint position in every time step. Fromthis point on, the postprocessing procedure is equivalent to experiments that track the motionof individual fluorescent tracer particles (e.g., [7,8,12]). Based on the discrete time sequenceof particle positions, the mean squared displacement (MSD) < ∆r p ( τ ) > is computed for anydesired time interval (that can be observed in a given simulation run) τ ∈ [ ∆t, t end ] as follows: < ∆r p ( τ ) > = 1 N τ N τ (cid:88) i =0 (cid:0) r p ( i∆t + τ ) − r p ( i∆t ) (cid:1) . (2.7)Here, N τ denotes the number of all distinct (but possibly overlapping) time intervals τ obtainedfor one simulation run. Given that the number of independent samples obtained for large timeintervals is naturally limited, only the first 10% of the maximal possible time intervals, i.e., τ ∈ [ ∆t, . · t end ] , will be considered in the statistical analyses. However, the remaining data isincluded 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 randomnetwork geometry as well as several realizations of the stochastic process of thermal forces will beconsidered. As the majority of the considered scenarios shows a subdiffusive behavior, the MSDcurves over the time interval will be presented and discussed instead of the (apparent) diffusionconstant, which obviously depends on the considered time interval and could still be computedfrom the MSD curves if desired. Other simulation results, such as the resulting axial strains ofthe fibers, will be presented and discussed for a few specific investigations wherever needed forinterpreting the system behavior.
3. 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 particlemobility. Second, the effect of additional attractive electrostatic interactions between the particleand the fiber network will be investigated in Section (b). Lastly, Section (c) analyzes the specialrole of fiber stiffness in the case of the most effective hindrance mechanism observed in thesimulations. (a) The effect of solely repulsive steric interactions
To begin with, the double-logarithmic plot of the particle’s MSD < ∆r p > as a function of thetime interval τ shown in Figure 4(a) confirms the validity of the applied Brownian dynamicsmodel, because the results obtained for free diffusion of the particle (red line with circles anderror bars) excellently match the analytical reference solution (black dashed line). Moreover, the r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... -3 -2 -1 -3 -2 -1
1 1 (a) (b)(c)
Figure 4.
Analysis of the particle mobility in presence of solely repulsive steric interactions with the fiber network. (a) Meansquared displacement (MSD) of the particle < ∆r p > as a function of the time interval τ : Mean and standard deviationover ten random realizations for the case of free diffusion (red line with circles and error bars) and three individual randomrealizations for a medium fiber volume fraction ¯ V f ≈ . (cyan lines). The analytical solution for the case of free diffusionis plotted as a reference (black dashed line). In addition, two sets (I) and (II) of three individual random realizations eachfor a very high fiber volume fraction ¯ V f ≈ . (green dotted lines) are included. To demonstrate the extreme effect ofcaged particles, all fibers have been fixed in space for all realizations with ¯ V f ≈ . . The gray background indicates therange of time intervals above 10% of the simulation time, where only few independent samples are available for computingthe MSD. (b) Network with very high fiber volume fraction ¯ V f ≈ . and overlay of all observed particle positions in onesimulation run corresponding to either the high MSD plateau value (I, blue) or the low value (II, orange). (c) Magnifieddetail showing the two compartments with irregular polygonal shape that the particle cannot leave. standard deviation of the mean over ten random realizations as indicated by the error bars isnegligible within the first 10% of the time interval range (indicated by the white background) thatis considered in the following analyses. Steric hindrance is insignificant as long as particle diameters are smaller than the smallestmesh sizes.
Turning to the influence of repulsive steric interactions, i.e., collisions between the particle andthe fiber network, Figure 4(a) shows that a medium fiber volume fraction ¯ V f ≈ . (cyan lines,see Figure 1(a) for an example of the network) has almost no perceptible influence on the MSD.Only above approximately τ = 0 . is a very subtle subdiffusive behavior observable for thethree individual realizations plotted here, which indicates the occurrence of a few collisionsif the particle travels over longer periods of time. On the one hand, this behavior of almostfree diffusion of the particle is expected from the range of cell diameters . − . µ m of this r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... irregular polygonal network compared to the particle diameter D p = 1 µ m . Moreover, previousexperimental studies have made very similar observations for almost neutral particles or highsalt concentrations that effectively shield any electrostatic interactions [7,8,12]. However, this isthe first computational study with a realistic, irregular polygonal fiber network geometry andthus an important confirmation that the effect of solely steric hindrance is indeed negligible inthis 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 compartmentsof random size and show confined diffusive behavior.
To push this to the limit where purely steric hindrance and thus filtering by particle sizewill become significant, a ten times larger number of Voronoi points n VP = 600 correspondingto a fiber volume fraction of ¯ V f ≈ . and cell diameters of . − . µ m , i.e., in the orderof the particle diameter D p = 1 µ m , have been applied. An example of the resulting networkarchitecture 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 tonormal diffusive behavior on very small time scales up to approximately τ = 0 .
05 s and eventuallyreaches a plateau value of < ∆r p > ≈ . µ m beyond τ ≈ . The other set (II), which are randomrealizations using the identical fiber network geometry, but a different initial position of theparticle, shows a significantly subdiffusive behavior already for the smallest considered timeintervals and quickly reaches a plateau value of < ∆r p > ≈ . µ m for any time interval longerthan τ ≈ . .This is an expected behavior for the diffusion of particles in an irregular network withcell diameters of the same order as the particle diameter and therefore randomly connectedsufficiently large cells that together form a polygonal volume surrounding the initial position ofthe particle, which the particle cannot leave under any circumstances. The specific compartmentsthat the particle is able to explore starting from either of the two initial positions in theidentical network are illustrated in Figure 4(b) and 4(c). Here, the overlay of all particlepositions throughout the entire simulation is shown in dark blue for the case of the higher MSDplateau value (I) and in orange for the case of the lower MSD plateau value (II) observed inFigure 4(a). This behavior is known as confined diffusion and has been theoretically described andexperimentally observed e.g., in the context of studying cadherin molecule mobility in plasmamembranes [55]. As discussed for instance in [1], such a filtering mechanism based on size clearlyhas a very effective selectivity and is applied by organisms to strictly preclude the access of anyobjects 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 inbiological hydrogels taking into account the continuous reorganization of biopolymer networksbased on both (de-)polymerization and the transient nature of crosslinks. Such a transportmechanism for relatively large particles has been observed in experiments [6] and recentlybeen investigated also in a theoretical and computational model [35], where crosslink bindingdynamics are influenced by the diffusing particle. It is also suggested that this kind of mechanismcould play a role in the selective permeability of the nuclear pore complex, for which thegoverning principles are still under debate (see e.g., the review article [1]). Replacing theVoronoi-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 futurestep. (b) The effect of additional attractive electrostatic interactions In addition to repulsive steric interactions, attractive electrostatic interactions due to uniformlydistributed, opposite charges on the particle and the fibers will now be considered. This has beenconfirmed to be the most effective hindrance mechanism for particles smaller than the mesh sizes r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... both in experiments (e.g., [7,8,12]) and simulations (e.g., [27,28]). Throughout this section, we thuskeep the fiber volume fraction fixed at the medium value ¯ V f ≈ . , which has been shown to nothave 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’ssurface charge Q p are shown in Figure 5 and compared to the analytical reference solution for freediffusion (black dashed line). As mentioned earlier, the fiber volume fraction ¯ V f ≈ . is identical -3 -2 -1 -3 -2 -1
1 1 (a) -3 -2 -1 -3 -2 -1
1 1 (b)
Figure 5.
Analysis of the hindrance of particle mobility due to attractive electrostatic interactions with the fiber network(in addition to repulsive steric interactions). Mean squared displacement (MSD) of the particle < ∆r p > as a function ofthe time interval τ : (a) Mean and standard deviation over five random network geometries and two random realizationseach for three different values of the particle’s surface charge: Low charge Q p = 0 . × − fC (green), mediumcharge Q p = 10 − fC (cyan), and high charge Q p = 8 × − fC (red). (b) All corresponding individual realizations.The medium fiber volume fraction ¯ V f ≈ . is identical for all of these realizations. The analytical solution for the caseof free diffusion is plotted as a reference (black dashed line). The gray background indicates the range of time intervalsabove 10% of the simulation time, where only few independent samples are available for computing the MSD. The bottomof the error bar is hidden for clarity wherever the corresponding value is negative. for all simulation runs. In particular, 5 different random network geometries with 2 randomrealizations of the stochastic forces each have been simulated for each of the three differentcharge values. Figure 5(b) shows these 10 independent realizations for each particle charge valueand the corresponding mean values and standard deviations are plotted in Figure 5(a). Onaverage, the degree of subdiffusion increases with the strength of attractive interaction, whichhas been suggested by previous experiments (e.g., [7,12]), and has been confirmed by previouscomputational studies using ordered (e.g., [26,27]) and unordered [28] arrays of straight, rigid,mutually orthogonal fibers. Here, the smallest charge value Q p = 0 . × − fC leads to a verysmall degree of subdiffusion, which is in fact quite similar to the one observed in the limit of nocharge shown in Figure 4(a). In contrast, medium (cyan) and high charge values (red), which area factor of 8 and 64 higher than the smallest charge value, significantly hinder the diffusion of theparticles already on (very) small time scales τ < . . In this regime, the slopes of the MSD curvesare however close to one, which suggests normal diffusive behavior with a decreased diffusionconstant. The variability of MSD values and slopes increases for longer time intervals, which indicatesthat particles randomly switch between distinct motion patterns.
This almost normal diffusive behavior with slope values close to one on very small time scaleschanges drastically for longer time intervals τ > . , where the individual realizations exhibit r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... slopes in the broad range from zero to one, and even some examples for superdiffusive behaviorwith a slope greater than one can be observed for medium charge values. Such a significantincrease in the variation of MSD values as well as slopes is in excellent agreement with theexperimental results from [12]. Their subsequent analysis of squared displacement values overtime revealed the existence of sudden trapping and escape events of a particle switching betweenan almost free, a loosely bound, and a tightly bound state. This hypothesis is confirmed by thesimulation results shown in Figure 6(a), where very similar motion patterns can be observed. (a)(b) (c) Figure 6.
Simulation results for one realization with high charge Q p = 8 × − fC and medium fiber volumefraction ¯ V f ≈ . that exhibits three distinct trapped states with sudden transitions (“jumps”) between them. (a) Squareddisplacement of the particle ∆r p over the simulation time t . (b) Particle trajectory with color (green, blue, pink, yellow,cyan and red) indicating the characteristic intervals of simulation time where the particle is either trapped (blue, yellowand red) or jumps between the trapped states (pink and cyan). (c) Detail view showing an overlay of the particle positionsin the three distinct trapped states (blue, yellow and red) and an overlay of the (oppositely charged) fiber aggregates thatare responsible for the effective trapping. The three colors (blue, yellow, and red) again match the corresponding timeintervals in (a) and (b). (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 ofopposite sign, has already been suggested in several experimental studies (e.g., [7,12]), the limitedspatio-temporal resolution of single particle tracking and imaging has so far precluded a directproof. The particle trajectory shown in Figure 6(b) and the detail view of the fiber aggregatescorresponding to the three observed trapped states illustrated in Figure 6(c) clearly confirmthis causality. Related observations have been made in the recently published computationalstudy [28], which for the first time considered disorder in the still rigid, straight, mutually r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... 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 biologicalhydrogels (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 inthe network that have a lateral spatial extension smaller than the particle diameter. Figure 6also suggests that the magnitude of the remaining thermal fluctuations observed in the squareddisplacement of the particle is a measure for the strength of the trapping and that this strengthis proportional to the local fiber, i.e., charge density. Finally, the results suggest that the strongerthe particle is immobilized, i.e., the smaller the fluctuations in the squared displacement are, thelonger the particle will remain at this location. The irregularity in the fiber/charge distribution gives rise to three distinct motion patterns, oneof which allows the particle to travel via successive jumps.
From a mechanical point of view, the spatially varying distribution of fibers and thus chargesgives rise to a random 3D potential field to be explored by the particle (similar to the one used inour 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 barrierand paths connecting them strongly depend on the characteristic geometrical properties of thenetwork such as connectivity, distribution of fiber segment lengths and the resulting distributionof cell and mesh sizes. Considering real biological hydrogels further extends the list of crucialinfluencing factors to the type and fractions of load-carrying components and their specificsurface charge distributions. All of these factors that shape the characteristic potential landscapewill 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 minimabeing restricted to a certain region, or (3) travel through the hydrogel via jumps – potentially alsoover long distances.All three motion patterns (1)–(3) can be identified among the realizations shown in Figure 5(b).For long time intervals τ > . , both patterns (1) and (2) lead to a similar behavior of confineddiffusion as obtained for the case of caged particles considered in the previous section . Theplateau MSD value thus allows to draw conclusions with respect to the volume enclosed by theparticle’s initial position and the location(s) of the potential minimum/minima (i.e., fiber/chargeagglomerations) that the particle has visited. By looking at the sequence of snapshots over theentire simulation time, it has been verified that those four (highly charged) particles with thesmallest MSD plateau values of < ∆r p > ≈ × − µ m are completely immobilized at onesingle location in the network (cf. motion pattern (1)). In contrast, those medium and highlycharged particles with intermediate values for the MSD slopes jump or smoothly transit betweenlocal fiber agglomerations in cycles (cf. motion pattern (2)) such as observed in the exampleshown in Figure 7. Finally, also motion pattern (3) has been identified in the set of individualrealizations, as already shown in Figure 6 for a particle with high surface charge. This case caneasily be identified among the MSD curves plotted in Figure 5(b) as the one with the largest MSDvalues for long time intervals.To conclude this section it can be stated that strong attractive forces mostly lead to motionpatterns (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 considerabledistances by means of a series of successive jumps. The effectiveness of this transport mechanismdepends 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 theparticle in the specific example shown in Figure 6 seems to be completely immobilized after twoconsecutive jumps, one might think of another special design of the fiber/charge distributionwith a certain degree of periodicity that could lead to an effective transport of the particle also Note however the difference in the diffusive behavior for small time intervals τ < . that allows to differentiate betweenpurely steric hindrance (cf. Figure 4(a)) and additional attractive interactions (cf. Figure 5(b)). r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... Figure 7.
Particle trajectory for one realization with medium charge Q p = 1 × − fC and medium fiber volumefraction ¯ V f ≈ . that smoothly cycles around a region of fiber/charge agglomeration. The continuous color schemeindicates the course of the simulation time t . over large distances and in small time intervals. Such a directed motion would maximize themobility of a diffusing particle.Altogether, the results of this section indicate that the irregularity of the potential landscapeexerts a crucial influence on the overall effectiveness of particle immobilization and thereforeon the selectivity of hydrogels. The aforementioned factors that shape the effective potentialfield explored by the particle are known to vary substantially between the multitude of differentbiological hydrogels (see e.g., [1]). Therefore, a systematic parametrization (and if required anextension) of the present computational model with respect to other classes of gels is expectedto be a valuable means for studying the species-specific variations of the general behavior andprinciples observed so far. In the long run, this might even lead to simulation-based predictiontools enabling a case-specific choice or design of drug delivery vehicles. Being able to optimizethe properties of the carrier to effectively attach the drug to the carrier and yet effectively diffusethrough the body would clearly be invaluable in this context. Brief discussion of computational aspects.
In total, 60 simulations with at least × time steps each have been conducted and evaluatedfor the results shown in this section (cf. Figure 5(b) for all realizations). One simulation typicallytook 2-4 days if run in parallel on 16 cores on a Linux cluster. The main drivers for thecomputational cost are the fiber volume fraction and the presence and strength of electrostaticinteractions leading to a complex interplay of repulsive steric and attractive electrostatic forcesthat require a fine temporal resolution (i.e., smaller time steps) and make the nonlinear problemmore 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 influencesthe amount of fiber deformation and – besides the realistic network geometry – is the secondunique feature of the present computational model. In the problem setup considered here,fiber deformations originate exclusively from contact and/or electrostatic interactions with theBrownian particle. Thus, the highest particle charge Q p = 8 × − fC will be considered in thissection, because the most frequent and strongest interactions can be expected in this case. Todetermine the point where fiber deformations begin to change the results, the value for Young’smodulus has been varied systematically starting from the theoretical limit of rigid fibers asoutlined already in Section (b). In this section, the results for rigid fibers will thus be compared to AMD Opteron 6128 r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... those obtained for a value of E = 0 . where the first differences can be observed, and to theones obtained for a ten times larger value E ∗ = 10 · E = 1 MPa , which are basically identical tothe case of rigid fibers. Note that – just as in the previous section – the fiber volume fraction willbe kept constant at the medium value ¯ V f ≈ . throughout this section.Figure 8 compares the MSD curves obtained for the low fiber stiffness resulting from Young’smodulus 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 -3 -2 -1 -3 -2 -1
1 1 (a) -3 -2 -1 -3 -2 -1
1 1 (b)
Figure 8.
Analysis of the influence of the fiber stiffness on the (hindrance of) particle mobility (due to attractive electrostaticand repulsive steric interactions with the fiber network). Mean squared displacement (MSD) of the particle < ∆r p > asa function of the time interval τ : (a) Mean and standard deviation over five random network geometries and two randomrealizations each for three different levels of the fiber stiffness: Low value for Young’s modulus E = 0 . (red, seealso Figure 5), medium value for Young’s modulus E ∗ = 1 MPa (yellow), and the limit of rigid fibers (cyan). (b) Allcorresponding individual realizations. The fiber volume fraction ¯ V f ≈ . and particle’s charge Q p = 8 × − fC are identical for all these realizations. The analytical solution for the case of free diffusion is plotted as a reference (blackdashed line). The gray background indicates the range of time intervals above 10% of the simulation time, where only fewindependent samples are available for computing the MSD. The bottom of the error bar is hidden for clarity wherever thecorresponding value is negative. fiber stiffness, the compared sets of ten realizations each are identical, in particular with respect tothe 5 different random network geometries and the 2 different sequences of the random stochasticforces each. The influence of fiber stiffness on particle mobility is insignificant within the range of stiffnessvalues reported for ECM gels.
There is no perceptible difference between the results for the medium fiber stiffness using E ∗ andthose for the limiting case of rigid fibers. Recall from Section (b) that the value E ∗ = 1 MPa isalready the lower limit of the wide range E = 1 MPa − of values reported for experimentswith collagen gels [48,49]. For more flexible fibers, we observe and speculate about a few competing effects that seem tocause an overall decrease of particle mobility.
First, let us look at the resulting overall effect in terms of the change in ensemble-averaged MSDcurves. For the lowest fiber stiffness E = 0 . considered in our simulations, the curves areshifted towards smaller MSD values (see Figure 8). The difference is not significant due to theconsiderable variability already observed in the results of the last section. However, these resultsallow to conclude that – for the given set of parameters – a value below E ∗ will begin to influencethe simulation results. As an explanation for this behavior, we speculate that the softer fibers r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... lead to an increased adhesive contact area because their shape adapts to the surface shape of theparticle and that this behavior in turn traps the oppositely charged particles more effectively. Inaddition, the softer fibers will presumably follow the particle’s thermal excitation more closelyand thus lead to smaller peak values of particle accelerations and interaction forces that actto separate the adhesive contact between fiber(s) and the particle, such that escape events willbecome less likely. Both effects might contribute to the observed overall reduction of particlemobility for fibers with a low stiffness. This reasoning is supported by the previously recognizedhigh importance of trapping and escape events as concluded from the results of the previousSection (b).Generally, there might be more mechanisms how very soft fibers can influence the particlemobility, including some that increase particle mobility. In an effort to shed some light onthis, we look at the direct, pair-wise comparison of a pair of simulations with different fiberstiffness values E and E ∗ and an otherwise identical setup, i.e., in particular, with identicalnetwork 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 -3 -2 -1 -3 -2 -1
1 1
Figure 9.
Detailed comparison between low fiber stiffness (red solid lines: E = 0 . ) and medium fiber stiffness(yellow dashed lines: E ∗ = 1 MPa ) for three individual random realizations (crosses, diamonds and circles). The plotshows the mean squared displacement (MSD) of the particle < ∆r p > as a function of the time interval τ . The analyticalsolution for the case of free diffusion is plotted as a reference (black dashed line). pair marked with crosses and second pair marked with diamonds) show nearly identical MSDcurves for E (red) and E ∗ (yellow), respectively. Conversely, the third pair of realizations withotherwise 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 cantrigger entirely different particle trajectories in individual realizations of the stochastic process.Given the already large variability in the results of individual random realizations for identicalfiber stiffness observed in the last section (see the right half of Figure 8(b)), this sensitivity is notsurprising. The stochastic nature of the problem impedes the direct investigation of differencesin particle behavior by means of a pair-wise comparison and underlines the importance of thealready 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 particlemotions as overlaying videos, we can still investigate a quite obvious mechanism that influencesthe particle mobility providing, in this case, nearly identical MSD curves. Indeed, one wouldintuitively argue that stiffer fibers constrain the motion of a particle more effectively as longas it remains in the state of adhesive contact with the fiber , because the softer fibers deform dueto 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 tworealizations. 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 r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... (a) (b) (c) Figure 10.
Comparison of the amount of fiber deformations visualized by an overlay of the configurations from all timesteps: (a) Medium fiber stiffness E ∗ preventing basically any deformations; (b) Ten times lower fiber stiffness E = 0 . E ∗ leading to noticeable, yet small fiber deformations. (c) Visualization of the magnitude and distribution of axial strains in thefiber network resulting from the thermal excitation of the particle for a low fiber stiffness E . time intervals τ ≥ . . However, these MSD curves also show that this effect is rather negligiblein terms of overall mobility of the particles. This is also supported by recognizing the smallmagnitude and localized extent of axial strains in the network (at an exemplarily chosen pointin time) in Figure 10(c).To conclude this first brief study of how deformations of semiflexible fibers influence themobility of oppositely charged particles, it can thus be stated that our results indicate a – rathercounter-intuitive – overall decrease in particle mobility if compared to (almost) rigid fibers. Wesuggest that this is the result of a decreased probability for escape events due to a) an increasein the adhesive contact area between particle and fiber and b) a higher ability of the fibers tofollow the thermal excitation of a trapped particle, thus reducing the peak values of particleaccelerations 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 rangeof stiffness values that have been reported for the collagen I fibers prevailing in ECM gels. Withinthis range, the fibers do not show any noticeable deformations for the scenario considered hereand the associated influence on the particle trajectory thus vanishes. Nevertheless, the influenceof fiber deformations is likely to be observed for other types of (biological) hydrogels with verythin or soft and therefore more flexible filaments such as mucin or F-actin, and maybe even forECM 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.
4. Conclusions and outlook
This article presents the first computational study of the diffusive mobility of particles inhydrogels with a realistic fiber network model. It proposes a novel computational approachbased 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 networkgeometries, and the beam-sphere interaction model for contact and electrostatic interactions.Following the validation of the model and the study of repulsive steric interactions only, theparticularly important effect of additional attractive electrostatic forces has been investigated.Finally, we have studied the role of fiber deformations in the latter case by means of additionalcomputational experiments. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... In the case of only repulsive steric interactions, it is found that the hindrance of the particlemobility is insignificant as long as the mesh sizes of the fiber network are larger than the particlediameter. If, however, the mesh sizes are in the order of the particle diameter, the particle iseffectively caged in a polygonal fiber hull of random shape and size and shows a behavior knownas 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 smallparticles, these are expected results validating the novel approach. However, this effect of sterichindrance will become highly relevant if the effective transport of relatively large particles asobserved in experiments [6] is considered. In this context, including the dynamics of the fibernetwork such as its self-assembly driven by the Brownian motion and transient reorganization ofcross-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 thenuclear 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 theoppositely charged particle, the prevailing notion that the degree of hindrance on averageincreases with the strength of attraction has been confirmed by the numerical experiments withfive different network geometries and two random realizations each. Moreover, an increasedvariability of the particle’s mean squared displacement values and slopes in the regime of longtime 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 prooffor the existence of distinct motion patterns of the particles, which explains the variability in theMSD curves. As hypothesized in the previous work [12], the particles stick to oppositely chargedfiber/charge aggregations experiencing more or less strong trapping and eventually escape due tothe 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 entire20s of simulation time, others smoothly or rapidly cycle between two local minima in the potentiallandscape. Both of these motion patterns lead to a behavior on longer time intervals that is verysimilar to the confined diffusion for caged, uncharged particles as described above. However, thediffusive mobility on short time intervals is significantly reduced as well due to sticking to thefibers. The third motion pattern observed is the one of several successive jumps that – at leasttheoretically – could serve as a transport mechanism also over longer distances if the potentiallandscape is formed accordingly and e.g., shows some degree of periodicity and directionalpreference.Altogether, these findings indicate that the precise shape of the effective 3D potential fieldexplored by the particle has a crucial influence on its mobility. In view of the broad varietyof biopolymer hydrogels with diverse chemical compositions and biophysical properties, thecurrent computational model could thus be leveraged to study the individual selective filteringbehavior for a large number of particle-hydrogel property combinations. Based on the recognizedimportance of the precise fiber/charge distribution in the system, two points seem to be ofparticular 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., byexperimentally analyzing the molecular architecture) and applied in the model. Second, thespecific composition and geometry of fiber networks should be determined (e.g., by processingelectron microscopy images) and applied in the model. Also the inclusion of the dynamic self-assembly and reorganization of networks mentioned above would be worth considering in thisrespect.As a last particular aspect investigated in this study, fiber deformations have been found tobe negligible within the range of realistic values for the stiffness of collagen I fibers prevailing inECM gels. To be more precise, varying the value for the Young’s modulus over the broad rangeof reported values for the considered ECM gels has led to identical results as obtained for thetheoretical limit of rigid fibers. If, however, more flexible fibers are considered, our simulationresults indicate an overall decrease of particle mobility if compared to (almost) rigid fibers – r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A .......................................................... an outcome that is rather counter-intuitive. We suggest that this is the result of a decreasedprobability for escape events due to a) an increase in the adhesive contact area between particleand fiber and b) a higher ability of the fibers to follow the thermal excitation of a trappedparticle, thus reducing the peak values of particle accelerations and interaction forces that actto 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 orF-actin, and maybe even for ECM gels as a result of dysregulated fiber stiffness. A more detailedinvestigation of this aspect, both in simulations and experiments, is thus considered a promisingavenue of future research.In addition to the presented simulation results and the gained insights, this study providesan extensive proof of concept for the application of the novel computational model. Asoutlined above, in the short to medium term many important findings especially for variousparticle/hydrogel-specific behaviors and mechanisms can be expected from applying differentparametrizations and, additionally, from integrating the suggested model extensions. In thelong term, further validation and advances of the present computational model toward a trulypredictive tool could ultimately lead to a case- and patient-specific choice or even design ofpharmaceuticals 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 isprovided or referenced in this article. The source code is part of a wider code base [41] that will be madeavailable in the future. The data is provided as electronic supplementary material.
Authors’ Contributions.
All authors contributed to writing the manuscript. WAW, OL and MJG designedthe research. MJG, CM and JFE developed the numerical model and methods. JK and MJG performed thenumerical experiments and analyzed the data.
Competing Interests.
All authors declare that there is no conflict of interest.
Acknowledgements.
We thank Kei W. Müller for useful discussions.
References
1. Witten J, Ribbeck K. 2017 The particle in the spider’s web: transport through biologicalhydrogels.
Nanoscale , 8080–8095.2. Bray J, Robinson GB, Byrne J. 1984 Influence of charge on filtration across renal basementmembrane films in vitro. Kidney International , 527–533.3. Dowd CJ, Cooney CL, Nugent MA. 1999 Heparan Sulfate Mediates bFGF Transport throughBasement Membrane by Diffusion with Rapid Reversible Binding. Journal of BiologicalChemistry , 5236–5244.4. Dellian M, Yuan F, Trubetskoy VS, Torchilin VP, Jain RK. 2000 Vascular permeability ina human tumour xenograft: molecular charge dependence.
British Journal of Cancer ,1513–1518.5. Olmsted SS, Padgett JL, Yudin AI, Whaley KJ, Moench TR, Cone RA. 2001 Diffusion ofMacromolecules and Virus-Like Particles in Human Cervical Mucus. Biophysical Journal ,1930–1937.6. Lai SK, O’Hanlon DE, Harrold S, Man ST, Wang YY, Cone R, Hanes J. 2007 Rapid transportof large polymeric nanoparticles in fresh undiluted human mucus. Proceedings of the NationalAcademy of Sciences , 1482 – 1487.7. Lieleg O, Baumgärtel RM, Bausch AR. 2009 Selective filtering of particles by the extracellularmatrix: an electrostatic bandpass.
Biophysical Journal , 1569–1577.8. Lieleg O, Vladescu I, Ribbeck K. 2010 Characterization of Particle Translocation throughMucin Hydrogels. Biophysical Journal , 1782–1789.9. Colwell LJ, Brenner MP, Ribbeck K. 2010 Charge as a Selection Criterion for Translocationthrough the Nuclear Pore Complex. PLOS Computational Biology , e1000747.10. Schuster BS, Suk JS, Woodworth GF, Hanes J. 2013 Nanoparticle diffusion in respiratorymucus from humans without lung disease. Biomaterials , 3439–3446. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
11. Xu Q, Boylan NJ, Suk JS, Wang YY, Nance EA, Yang JC, McDonnell PJ, Cone RA, Duh EJ,Hanes J. 2013 Nanoparticle diffusion in, and microrheology of, the bovine vitreous ex vivo.
Journal of Controlled Release , 76–84.12. Arends F, Baumgärtel R, Lieleg O. 2013 Ion-Specific Effects Modulate the Diffusive Mobilityof Colloids in an Extracellular Matrix Gel.
Langmuir , 15965–15973.13. Zhang X, Hansing J, Netz RR, DeRouchey JE. 2015 Particle Transport through Hydrogels IsCharge Asymmetric. Biophysical Journal , 530–539.14. Käsdorf BT, Arends F, Lieleg O. 2015 Diffusion Regulation in the Vitreous Humor.
BiophysicalJournal , 2171–2181.15. Abdulkarim M, Agulló N, Cattoz B, Griffiths P, Bernkop-Schnürch A, Borros SG, GumbletonM. 2015 Nanoparticle diffusion within intestinal mucus: Three-dimensional response analysisdissecting the impact of particle surface charge, size and heterogeneity across polyelectrolyte,pegylated and viral particles.
European Journal of Pharmaceutics and Biopharmaceutics , 230–238.16. Arends F, Chaudhary H, Janmey P, Claessens MMAE, Lieleg O. 2017 Lipid Head GroupCharge and Fatty Acid Configuration Dictate Liposome Mobility in Neurofilament Networks. Macromolecular Bioscience , 1600229.17. Johansson L, Löfroth JE. 1993 Diffusion and interaction in gels and solutions. 4. Hard sphereBrownian dynamics simulations. The Journal of Chemical Physics , 7471–7479.18. Saxton M. 1994 Anomalous diffusion due to obstacles: a Monte Carlo study. Biophysical Journal , 394–401.19. Netz PA, Dorfmüller T. 1997 Computer simulation studies of diffusion in gels: Modelstructures. The Journal of Chemical Physics , 9221–9233.20. Pei H, Allison S, Haynes BMH, Augustin D. 2009 Brownian Dynamics Simulation of theDiffusion of Rods and Wormlike Chains in a Gel Modeled as a Cubic Lattice: Applicationto DNA.
The Journal of Physical Chemistry B , 2564–2571.21. Stylianopoulos T, Diop-Frimpong B, Munn LL, Jain RK. 2010 Diffusion Anisotropy inCollagen Gels and Tumors: The Effect of Fiber Network Orientation.
Biophysical Journal ,3119–3128.22. Kamerlin N, Elvingson C. 2016 Tracer diffusion in a polymer gel: Simulations of static anddynamic 3D networks using spherical boundary conditions. Journal of Physics CondensedMatter .23. Saxton M. 1996 Anomalous diffusion due to binding: a Monte Carlo study. Biophysical Journal , 1250–1262.24. Zhou H, Chen SB. 2009 Brownian dynamics simulation of tracer diffusion in a cross-linkednetwork. Physical Review E , 21801.25. Stylianopoulos T, Poh MZ, Insin N, Bawendi MG, Fukumura D, Munn LLL, Jain RK.2010 Diffusion of Particles in the Extracellular Matrix: The Effect of Repulsive ElectrostaticInteractions. Biophysical Journal , 1342–1349.26. Miyata T. 2012 Brownian Dynamics Simulation of Self-Diffusion of Ionic Large SoluteMolecule in Modeled Polyelectrolyte Gel. Journal of the Physical Society of Japan , SA010.27. Hansing J, Ciemer C, Kim WK, Zhang X, DeRouchey JE, Netz RR. 2016 Nanoparticle filteringin charged hydrogels: Effects of particle size, charge asymmetry and salt concentration. European Physical Journal E , 53.28. Hansing J, Netz RR. 2018 Particle Trapping Mechanisms Are Different in Spatially Orderedand Disordered Interacting Gels. Biophysical Journal , 2653–2664.29. Humphries DL, Grogan JA, Gaffney EA. 2017 Mechanical Cell-Cell Communication in FibrousNetworks: The Importance of Network Geometry.
Bulletin of Mathematical Biology , 498–524.30. Cyron CJ, Wall WA. 2010 Consistent finite-element approach to Brownian polymer dynamicswith anisotropic friction. Physical Review E , 66705.31. Slepukhin VM, Grill MJ, Müller KW, Wall WA, Levine AJ. 2019 Conformation of a semiflexiblefilament in a quenched random potential. Physical Review E , 042501.32. Cyron CJ, Müller KW, Schmoller KM, Bausch AR, Wall WA, Bruinsma RF. 2013a Equilibriumphase diagram of semi-flexible polymer networks with linkers. Europhysics Letters , 38003.33. Cyron CJ, Müller KW, Bausch AR, Wall WA. 2013b Micromechanical simulations ofbiopolymer networks with finite elements.
Journal of Computational Physics , 236–251.34. Müller KW, Bruinsma RF, Lieleg O, Bausch AR, Wall WA, Levine AJ. 2014 Rheology ofSemiflexible Bundle Networks with Transient Linkers.
Physical Review Letters , 238102. r s pa . r o y a l s o c i e t y pub li s h i ng . o r g P r o c R S o c A ..........................................................
35. Goodrich CP, Brenner MP, Ribbeck K. 2018 Enhanced diffusion by binding to the crosslinks ofa polymer gel.
Nature Communications , 4348.36. Marczynski M, Käsdorf BT, Altaner B, Wenzler A, Gerland U, Lieleg O. 2018 Transient bindingpromotes molecule penetration into mucin hydrogels by enhancing molecular partitioning. Biomaterials Science , 3373–3387.37. Meier C, Popp A, Wall WA. 2019 Geometrically Exact Finite Element Formulations for SlenderBeams: Kirchhoff-Love Theory Versus Simo-Reissner Theory. Archives of ComputationalMethods in Engineering , 163–243.38. Meier C, Grill MJ, Wall WA, Popp A. 2018 Geometrically exact beam elements and smoothcontact schemes for the modeling of fiber-based materials and structures. International Journalof Solids and Structures , 124–146.39. Grill MJ, Wall WA, Meier C. 2020 A computational model for molecular interactions betweencurved slender fibers undergoing large 3D deformations with a focus on electrostatic, van derWaals, and repulsive steric forces.
International Journal for Numerical Methods in Engineering , 2285–2330.40. Arends F, Nowald C, Pflieger K, Boettcher K, Zahler S, Lieleg O. 2015 The biophysicalproperties of basal lamina gels depend on the biochemical composition of the gel.
PLoS ONE , e0118090.41. BACI: A Comprehensive Multi-Physics Simulation Framework. 2020https://baci.pages.gitlab.lrz.de/website. .42. Rycroft C. 2009 Voro++: a three-dimensional Voronoi cell library in C++. Technical reportUnited States.43. Blender Foundation. 2019 Blender 2.80. .44. Yurchenco PD, Ruben GC. 1987 Basement membrane structure in situ: evidence for lateralassociations in the type IV collagen network.. The Journal of Cell Biology , 2559–2568.45. Reissner E. 1981 On finite deformations of space-curved beams.
Zeitschrift für AngewandteMathematik und Physik (ZAMP) , 734–744.46. Simo JC. 1985 A Finite Strain Beam Formulation. The Three-Dimensional Dynamic Problem.Part I. Computer Methods in Applied Mechanics and Engineering , 55–70.47. Simo JC, Vu Quoc L. 1986 A Three Dimensional Finite Strain Rod Model Part II:Computational Aspects. Computer Methods in Applied Mechanics and Engineering , 79–116.48. Jansen KA, Licup AJ, Sharma A, Rens R, MacKintosh FC, Koenderink GH. 2018 The Role ofNetwork Architecture in Collagen Mechanics. Biophysical Journal , 2665–2678.49. van der Rijt JAJ, van der Werf KO, Bennink ML, Dijkstra PJ, Feijen J. 2006 MicromechanicalTesting of Individual Collagen Fibrils.
Macromolecular Bioscience , 697–702.50. Jeleni´c G, Crisfield MA. 1999 Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics. Computer Methods in Applied Mechanics andEngineering , 141–171.51. Crisfield MA, Jeleni´c G. 1999 Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation.
Proceedings of the RoyalSociety of London. Series A: Mathematical, Physical and Engineering Sciences , 1125–1147.52. Grill MJ, Meier C, Wall WA. 2019 Investigation of the peeling and pull-off behavior of adhesiveelastic fibers via a novel computational beam interaction model.
The Journal of Adhesion,https://doi.org/10.1080/00218464.2019.1699795 .53. Doi M, Edwards SF. 1988
The theory of polymer dynamics . Oxford university press.54. Meier C, Popp A, Wall WA. 2016 A finite element approach for the line-to-line contactinteraction of thin beams with arbitrary orientation.
Computer Methods in Applied Mechanicsand Engineering , 377–413.55. Kusumi A, Sako Y, Yamamoto M. 1993 Confined lateral diffusion of membrane receptorsas studied by single particle tracking (nanovid microscopy). Effects of calcium-induceddifferentiation in cultured epithelial cells.
Biophysical Journal65