Magdalini Ntetsika and Panayiotis Papadopoulos

Department of Mechanical Engineering,University of California,Berkeley,CA,94720-1740,USA

ABSTRACT A new and computationally efficient version of the immersed boundary method,which is combined with the coarse-graining method,is introduced for modeling inextensible filaments immersed in low-Reynolds number flows.This is used to represent actin biopolymers,which are constituent elements of the cytoskeleton,a complex network-like structure that plays a fundamental role in shape morphology.An extension of the traditional immersed boundary method to include a stochastic stress tensor is also proposed in order to model the thermal fluctuations in the fluid at smaller scales.By way of validation,the response of a single,massless,inextensible semiflexible filament immersed in a thermally fluctuating fluid is obtained using the suggested numerical scheme and the resulting time-averaged contraction of the filament is compared to the theoretical value obtained from the worm-like chain model.

KEYWORDS Semiflexible biopolymers;immersed boundary method;coarse-graining;actin filaments;fluid-structure interaction;thermal fluctuations;persistence length

1 Introduction

Living cells display a high degree of internal mechanical and functional organization and their intracellular biopolymeric scaffold,the cytoskeleton,plays a key role in that [1].The actin cortex is the part of the cytoskeleton attached to the cell membrane in most eukaryotic cells.The cortical actin cytoskeleton plays a fundamental role in cell shape,which is maintained through structural stiffness and rheology.Actin cortex is a thin actomyosin network that underlies the plasma membrane,consisting of actin filaments cross-linked by actin-binding proteins and containing motor proteins that generate stress within the network [2].In order for the actin filaments to become cross-linked,they undergo thermal fluctuations to find a cross-linking partner and when cross-links between actin filaments are formed,the amplitudes of filament fluctuations are reduced,supporting the formation of additional cross-links [3].

Actin filaments are biopolymers with sufficient contour length to exhibit significant thermal bending fluctuations,in the order of approximately 1% of their contour length.However,their diameter can be as large as ten nanometers or more,giving them noteworthy bending rigidity.Thus,actin filaments are said to be semiflexible in the sense that their bending stiffness is large enough for the bending energetics—which favors a straight conformation—to just out-compete the entropic tendency of a chain to crumple up into a random coil [4].Therefore,semiflexible polymers exhibit small,yet significant,thermal fluctuations around a straight conformation.Furthermore,the semiflexible filaments are practically inextensible,i.e.,their backbone cannot be stretched or compressed.On the scale of several nanometers to micrometers,biopolymers are often effectively modeled as inextensible elastic rods or fibers with finite resistance to bending [5].This is the essence of the classical worm-like chain (WLC) model by Kratky et al.[6].The competition between entropic and energetic effects in semiflexible polymers gives rise to many physical properties and the semiflexible nature of the actin filaments also has major implications on how they interact with each other to form cross-linked networks [7,8].One may think of a single actin filament as a chain that can respond to forces or thermal fluctuations by bending and end-to-end compression relative to its full contour length.

There has been deep interest in studying the mechanical response of biological tissues the past decades,and more specifically,in understanding the mechanical properties of biopolymeric networks,since they play an important role in cell motility [9,10] and mechanotransduction[11–13].Furthermore,the investigation of the behavior of semiflexible filaments in viscous shear flows at low Reynolds numbers has also gained a lot of interest due to the relevance of its applications in areas involving biological systems like DNA [14,15],polymers [16] and proteins,but also in areas such as biotechnology that involve natural and synthetic fibers [17].These fluidstructure interaction (FSI) problems are quite challenging because of the complex interplay of hydrodynamic stresses and the corresponding fiber conformity,however,a two-way interaction between the immersed filaments and the surrounding fluid is tantamount in order to get a better understanding of the underlying physical behavior.

Yamamoto et al.[18] proposed a method for simulating the dynamic behavior of rigid and flexible fibers in a flow field,with the fibers regarded as made up of spheres that are lined up and bonded to each neighbor.Tornberg et al.[19] studied the dynamics of slender filaments suspended in Stokesian fluids employing a non-local slender body theory.In most of these studies,the hydrodynamic interaction was neglected,therefore,there was limited information about the underlying FSI.The fluid was considered as a passive medium and coupling between fluid and structure was one-way.The Immersed Boundary Method (IBM) of Peskin [20],which accounts for the two-way interaction between filament and fluid,has gained substantial popularity the past years in studying the dynamics of filaments in flow fields.However,in most of the studies where the dynamics of inextensible semiflexible filaments were studied using IBM,as in [21,22],a very high stretching stiffness was used for the filaments to approximate inextensibility,which significantly restricted the time step,thus increasing substantially the computational cost.Wiens et al.[23] used a generalized IB method which can be viewed as a type of penalty method in which the rod is only approximately inextensible.Similarly,Huang et al.[24] used a modified IBM with an extra inextensibility condition to strictly enforce the filament inextensibility.In addition,Kim et al.[25] introduced a penalty immersed boundary method to simulate the dynamics of inextensible vesicles in an incompressible viscous fluid by using two Lagrangian immersed boundaries connected through stiff springs to represent the real immersed boundary for different purposes.Also,Ong et al.[26] developed an immersed boundary projection method based on an unconditionally energy stable scheme to simulate the vesicle dynamics in a viscous fluid.

In the preceding studies the inextensibility constraint is not enforced strongly and this may lead to numerical errors [24] that cause the filament length to vary over time.On the other hand,when inextensibility is enforced by Lagrange multipliers,an extended nonlinear system of equations needs to be solved at every time step,thus again increasing substantially the computational time.An alternative primal approach,more straightforward and faster computationally,to deal with the inextensibility of semiflexible filaments is to use the Coarse-Graining Method(CGM) by Moreau et al.[27].CGM uses discrete models,where the filament is partitioned into a discrete number of straight segments and the elastic interaction coupling neighboring nodes/joints is described via discrete elastic connectors encoding the filament’s resistance to bending.However,up to now,the CGM for filaments has been used by neglecting the two-way interaction between the filament and the fluid.

In this study,a new computationally efficient version of the IBM,which is combined with the CGM,is introduced in this study for modeling inextensible filaments in low-Reynolds number flows.An extension of the traditional IBM to include a stochastic stress tensor is also proposed in order to model the thermal fluctuations in the fluid in smaller scales.The proposed numerical scheme is validated by comparing the response of a single actin filament immersed in a thermally fluctuating fluid to the theoretical values obtained using the WLC model.

The remainder of the article is organized as follows.The WLC model is reviewed in Section 2 and the theoretical value of the time-averaged contraction for a single inextensible filament under hydrodynamic thermal fluctuations is derived.The mathematical formulation of the coupled system filament-fluid and the suggested numerical procedure are described in Section 3.In Section 4,the behavior of a single,massless,inextensible,and semiflexible filament immersed in a thermally fluctuating fluid is obtained using the suggested numerical scheme and the resulting time-averaged contraction of the filament is compared to the theoretical value obtained from the WLC model for the sake of validation.This is followed by a concluding reflection on the findings in Section 5.

2 Theoretical Background:The Worm-Like Chain Model

The mechanical behavior of semiflexible filaments is usually described by the WLC model [6].This assumes that the filament is inextensible,has linearly elastic bending energy,and is subjected to thermal fluctuations.A homogeneous,inextensible and semiflexible filament of straight lengthLand circular cross section of radiusais taken to be fixed at one end and free at the other.The domain of the filament is parametrized by its arc-lengthsand the position of a typical point is denoted y(s).The motion of the filament is assumed to be confined to a plane and its bending energy is given by

whereκ=is the (linearized) curvature andκb=EIis the flexural rigidity,expressed in terms of the Young’s modulusEand the second moment of areaI=Using the Equipartition Theorem [28,Chapter 7] one may calculate the thermal average angular correlation between distant points along the filament,for which

where t(s)=is the (unit) tangent vector along the filament,kBis the Boltzmann constant,andTis the absolute temperature of the fluid.The persistence lengthlpis a characteristic length over which the relative orientation of two such tangent vectors along the filament remains correlated,and is defined asfor a planar motion [1,Section 3.2.3] [5].In simple terms,a filament in thermal equilibrium inside a fluid will appear rather straight over lengths that are short compared with this persistence length,while it will begin to exhibit a random,contorted shape otherwise.

For the filament under consideration here,it is assumed thatL≤lp,therefore it is expected to remain nearly straight with small transverse fluctuations.Let thex-axis define the average orientation of the filament segment,and leturepresent the transverse displacement taken to be a function ofxand timet.The functionu(x,t) may be represented by a Fourier series as

whereqare the wave numbers defined asq=nπ/L,wheren=1,2,3,...,anduqare the corresponding amplitudes.Such a representation is appropriate for the case of a nearly straight filament with boundary conditionu=0 atx=0 and no restraint atx=L,as in Fig.1.Since the transverse displacement is assumed infinitesimal,the local orientation and curvature of the filament are given byand,respectively.In view of (1) and (3),the bending energy is given by

as also in [29].

Figure 1:Configuration of filament fixed at one end,where u(x;t) denotes the transverse displacement u from an initial straight line (dashed).The magnitude of u is exaggerated for illustration

Since the filament is inextensible,the total arc-length of the filament remains unchanged under the influence of the fluctuations.Thus,the arc length dsof a short segment is approximately given byThe contraction of the chain relative to its full contour length in the presence of thermal fluctuations inuis then

where use is made of the approximationfor anyα,as well as of the absence of distinction between integrals over the arc-length and its projection on thex-axis.In view of (3),the contraction in (5)3is expressed as

By way of background,recall that the ensemble average＜A(t)＞of an observable quantityA,which is a function of a variableΓat timet,is defined as

see,e.g.,[30,Chapter 3].Here,γis the phase space ofΓandf(Γ,t) is the phase-space distribution function,defined as the probability that the system will be in stateΓat timet.The probability density function satisfies the standard consistency condition

In the context of the present problem,if the filament is in equilibrium at temperatureT,the ensemble average value＜u2q＞of the amplitude-squared is then deduced in accordance with the canonical distribution [28,Chapter 6],and is expressed as

where “ln” denotes the natural logarithm.A straightforward calculation results in

therefore,for anyq,the ensemble average of the amplitude-squared is written as

in terms of the original lengthLand the persistence lengthlp.Likewise,using (6) and (11),the ensemble average＜Δ L＞of the contractionΔ Lis expressed as

It can be observed from (11) that the ensemble average of the amplitude-squared of the bending fluctuations diminishes rapidly for higher-order modes due to the (negative) fourth-power dependence on the wave number.Also,(12) implies that the longer the persistence length,the smaller the ensemble average of the contractionΔL.

The probability density functionfin (7) is independent rendering the ensemble stationary.In this case,the ergodic hypothesis [28,Chapter 15] states then that the ensemble average over all accessible systems is equal to the time average over a large number of observation of a single system.That is,given any observable quantityAin a stationary ensemble,

for sufficiently large observation timet.Therefore,the ensemble averages＜u2q＞and＜ΔL＞in Eqs.(11) and (12),respectively,can be estimated by such time averaging.

3 Mathematical Formulation and Numerical Procedure

The fluid and the immersed semiflexible filament constitute a coupled mechanical system.The inextensible filament’s motion is driven by the fluid’s velocity field,while,at the same time,the filament exerts force on the fluid,thus affecting its motion.The equations of motion that describe the coupled system are derived and discussed in the remainder of this section.

3.1 The Asymptotic Coarse-Grained Elastohydrodynamics

Consider an inextensible massless filament of lengthLimmersed in an incompressible viscous fluid,and recall that the position of a point of the rod be denoted by y(s).The filament is embedded in a two-dimensional space associated with orthonormal basis {ex,ey} and is subjected to external contact force fh(s) per unit length due to the hydrodynamic interactions.

Assuming quasi-static loading conditions,the equilibrium equations are written as

where n and m are the (internal) axial force and moment sustained by the filament.Integrating the force equilibrium Eq.(14)1over the entire filament leads to

where it is further assumed that the boundary forces vanish,that is,n(0)=n(L)=0.

Following Moreau et al.[27],the filament is partitioned intoNstraight subdomains (elements)of lengthΔs=L/N.The Frenet basis along the filament is given by the unit tangent and normal vectors (e‖,e⊥),respectively,as shown in Fig.2.Also,the angle between e||and exis denotedθ.The discrete counterpart of the equilibrium Eq.(15) takes the form

where Fhirepresents the resultant external force experienced by thei-th element.For a filament free of boundary moments,that is,assuming,m(L)=m(0)=0,one may take the integral of(14)2over the entire filament,use integration by parts,and invoke (14)1and (16) to conclude that

where Mi,y0is the moment of the external force acting on thei-th element about point y0=y(0).Upon integrating Eq.(14)2over the domain ((j-1)Δs,L),j=2,...,N,and taking into account(17),it follows that

where mj=m((j-1)L/N),j=2,...,N,is the moment at the left end-point of thej-th element.Given that the moment at any point it is defined aswhereis the curvature and ez=ex×ey,one may use the backward finite difference formula to find that

whereαj=θj-θj-1,j=1,2,...,N,is the angle between ei-1,‖and ei,‖,as shown in Fig.2,and thus

withα0=θ0.

Figure 2:Parametrization of the continuous and discretized filament

In the low-Reynolds number regime,the hydrodynamic force experienced by the filament immersed in fluid with velocity field v can be defined according to the Resistive Force Theory [31]as

where v(y) is the fluid velocity interpolated at the filament positions,andξandηare the normal and tangential drag coefficients,respectively.Here,the external contact force fh(s) is due to the resistance that the filament points experience as they move in the fluid.

The position y at timetis a point in the filament that satisfies y(s)=y0+e‖(τ)dτ,so that in the discrete case

fori=1,...,N,thus yi=y((i-1)L/N) thus satisfying the inextensibility constraint from the outset.Here,ek,‖denotes the tangent Frenet vector in thek-th element.

There are nowN+2 parameters describing the position of the filament,that is,(y0,α0,...,αN-1).To determine them,there are two total force balance equations in (16),one torque balance in (17),andN-1 equations for the internal moment balance in (18),thus rendering this elastohydrodynamic system closed.

With slight abuse of notation,let yi(s) denote the current position of a filament point on thei-th element,so that

The total hydrodynamic force on thei-th element is found from (21),with the aid of (16) and(23),to be

Note that,in general,v varies along the domain (iΔs,(i+1)Δs).However,to within a small error,the velocity is assumed constant in each element.Using Eqs.(17) and (21),one finds that

Therefore,Eqs.(16)–(18),with the aid of (24) and (25),take the form

These coarse-grained elastohydrodynamics equations,in conjunction with Eqs.(20) and (22),can be cast a system of ordinary differential equations of the form

where [Y]=[x0y0θ0α1...αN-1]Tare theN+2 parameters with (x0,y0)=y0,and [V] is a vector of size 3N,containing in the first 2Nentries thex-andy-components of the flow velocity at the filament pointsi=0,1,...N-1 with the remaining entries being zero.

The matrix [A] has dimensions (N+2)×3Nand its coefficients are given,for alli,j=0,...,N-1,by Ifj＜i,thenai+2,j=ai+2,N+j=ai+2,2N+j=0.Also,the column vector [B] of sizeN+2 is given by

In addition,matrix [Q] is a 3N×(N+2) transformation matrix defined as

where [Q1] and [Q2] areN×Nmatrices whose elements are given by the general formula

The system of equations in (27) is integrated in time using an implicit first-order method with constant time stepΔt.Writing the system of equations at timetn+1,one can obtain [Yn+1] by solving implicitly the system of equations

in the typical time domain (tn,tn+1].

3.2 Hydrodynamic Fluctuations and Equations of Motion

As one approaches smaller length scales,in the order ofμm,thermal fluctuations play an essential role in the description of the fluid flow.Thermal fluctuations may be included in the continuum description of the fluid by means of additional stochastic fluxes.The resulting equations of motion for the fluctuating fluid turn out to be stochastic partial differential equations.Landau et al.[32,Chapter 9] proposed such equations,which include an additional stochastic stress tensor in the Navier-Stokes equations,the so-called Landau-Lifshitz Navier-Stokes (LLNS)equations.

To account for thermal fluctuations,the Cauchy stress tensor T for an incompressible viscous fluid can be modified as

wherepis the pressure,μis the dynamic viscosity of the fluid,D is the rate-of-deformation tensor,andstands for the stochastic stress tensor,which models the inherent molecular fluctuations in the fluid.The required stochastic properties ofhave been derived by Landau et al.[32] in order to satisfy the fluctuation dissipation theorem [32],according to

whereδ(·) is the Dirac delta function,I is the symmetric fourth-order identity tensor,＜·＞denotes ensemble average and ⊗denotes tensor product.

The low-Reynolds number Navier-Stokes equations for an incompressible,Newtonian fluid with the additional stochastic stress tensor to account for the thermal fluctuations can be written,in the absence of body forces,as

where,again,v is the fluid’s velocity field and alsoρis the density of the incompressible fluid.As it turns out,the flow in cytosol is of low Reynolds number,in the order of 10-9,in large part due to the cell’s size,which is in the order ofμm [1].For such flows it is reasonable to adopt the unsteady Stokes approximation,where the convective rate of changein the acceleration is neglected,but the spatial time derivative termis retained in order to capture the effects of the unsteady flow.

Consider now an incompressible viscous fluid occupying a two-dimensional domainΩand undergoing thermal fluctuations.The IBM formulation in its strong form can be understood as an enrichment of the two-dimensional Navier-Stokes equations accounting also for the forces generated by the deformation of the immersed body,with the linear momentum balance equations in (35)2for the fluid taking the form

where F is the force that the filament exerts to the fluid.

In the discrete case,the force term F included in the Navier-Stokes equations is defined to be equal and opposite to the equivalent hydrodynamic force Fhiexperienced by the fiber moving with the fluid with velocity field according to the Resistive Force Theory [31] and Eqs.(21) and (24).Also,the perpendicular and parallel drag coefficients in Eq.(21) for a filament of cross-sectional radiusrand length 2lare defined,following Lighthill’s classical analysis [33],as

There exist many choices for representing the stochastic stress tensorin Eq.(36)1[34].As suggested in [35–38],a formulation for the stochastic stress tensorthat requires the fewest possible random numbers and satisfies the fluctuation dissipation theorem [32] and Eq.(34),is expressed as

whereΔV==vol(Ωe),whereVis the volume of the fluid domainΩ,nelis the total number of elements in the Eulerian spatial discretization of the fluid,andΩeis the domain of elementeused in this (uniform) discretization.The tensoris symmetric,defined aswhere a realization of R is sampled using a stream of independent,standard normally distributed random numbers at each time stepΔtfor every elementΩeof the fluid domain.

3.3 Summary of the Numerical Algorithm for the Hybrid IBM-CGM

The implementation of the proposed numerical algorithm for simulating flexible inextensible filaments immersed in a fluid can be summarized as follows:

Algorithm 1 Hybrid IBM-CGM algorithm Initialize fluid field variables v0 and [p0]Initialize immersed curve position [y0]Set n=0 While n＜n_step do Compute [Vn] by interpolating [vn] at the filament positions Compute the filament position [yn+1] using Eq.(32)Compute [Fnhi]=[An][Vn]Distribute Lagrangian interaction force [Fn] to the Eulerian grid Calculate [images/BZ_178_358_2092_386_2148.png] using Eq.(38) for every cell of the fluid Solve for [vn+1] and [pn+1] based on Eq.(36)n←n+1 end while

4 Model Validation

This section focuses on the mechanical response of a single massless filament immersed in a thermally fluctuating fluid and aims to validate the results of the computational model described in the previous section by comparison with the theoretical prediction based on the WLC model in Section 2.

The biopolymers that comprise the cytoskeleton consist of aggregates of large globular proteins that are bound together rather weakly,as compared with most synthetic,covalently bonded polymers [5].Nevertheless,they can be surprisingly strong,due to their relatively large diameter,which makes their bending rigidity the dominant attribute determining their mechanical behavior on the cellular scale.Even with this mechanical resistance to bending,however,cytoskeletal filaments can still exhibit significant thermally induced bending fluctuations because of Brownian motion in the surrounding fluid.Actin filaments,owing to their relatively large bending stiffness,have long persistence length compared to their total contour length,which practically means that an actin filament in thermal equilibrium in a fluid will appear rather straight.Here,the results of simulations of an inextensible filament immersed in a fluid domain with thermal noise are used to verify that,after a sufficient amount of time,the time-averaged contraction＜ΔL＞tof the filament approaches the theoretical value of the ensemble average＜ΔL＞derived in Section 2.

The inextensible filament is modeled by means of the CGM described in Section 3.1.The filament has contour length ofL=5μm (typical lengths range from 100 nm to a few microns [39])and is initially straight and fixed at one end,as shown in Fig.1.Its bending rigidity is taken equal toκb=10-19Nm.The filament is immersed in a fluid domain which has temperature ofT=300 K,dynamic viscosityμ=0.001 Pa·s-1,and fluid densityρ=1000 kg/m3,since cytosol is composed mainly of water.The fluid domain is of size 10-5μm×10-5μm and the stochastic Navier-Stokes equations in (36) are solved in it under periodic boundary conditions.The fluid domain is discretized using 40×40 Q2-Q1 Taylor-Hood elements [40],while the filament is discretized byN=10 elements.For both the solid and the fluid,a constant time stepΔt=10-5s is used.The algorithm described in Section 3.3 is employed to resolve the fluid-solid interaction.

Fig.3 depicts the results of the simulation described above.The time-averaged contraction＜ΔL＞tis plotted as a function of the number of time steps.It can be concluded from the plot that after around 6 * 106time steps,＜ΔL＞tconverges to a value which is approximately equal to 3.31×10-13.This compares very well to the theoretical value for the ensemble average obtained from Eq.(12),which is approximately equal to 3.45×10-13for a relative error of approximately 4%.In view of this agreement,the comparison may serve as validation of the proposed numerical algorithm used to simulate the fluid-structure interaction of the immersed filaments under thermal fluctuations.

Figure 3:Time average (log scale) of inextensible filament’s contraction＜ΔL＞t as a function of the number of time steps under thermal fluctuations

5 Conclusions

In this study,a modified and computationally efficient version of the Immersed Boundary Method,combined with the Coarse-Graining Method,was proposed for modeling inextensible semiflexible filaments in low-Reynolds number flows.Thermal fluctuations in the fluid were modeled by including a stochastic stress.The mechanical behavior of a massless,inextensible,and semiflexible filament immersed in a thermally fluctuating fluid was investigated using the suggested method.The resulting time-averaged contraction of the filament compares very favorably to the theoretical value for the ensemble average of the same quantity,as obtained from the Worm-like Chain model.On the basis of this analysis,the proposed hybrid algorithm appears to be both robust and accurate,and could offer a reliable means for investigating the combined effect of multiple (and possibly interacting) filaments in low-Reynolds number flows.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.