Segmental Lennard-Jones Interactions for Semi-flexible Polymer Networks
Carlos Floyd, Aravind Chandrasekaran, Haoran Ni, Qin Ni, Garegin A. Papoian
SSegmental Lennard-Jones Interactions for Semi-flexible PolymerNetworks
Carlos Floyd , Aravind Chandresekaran , Haoran Ni , Qin Ni , and Garegin A.Papoian Biophysics Program, University of Maryland, College Park, MD 20742 USA Department of Chemistry and Biochemistry, University of Maryland, College Park,MD 20742 USA Department of Chemical and Biomolecular Engineering, University of Maryland,College Park, MD 20742 USA Institute for Physical Science and Technology, University of Maryland, College Park,MD 20742 USAFebruary 24, 2021
Abstract
Simulating soft matter systems such as the cytoskeleton can enable deep understanding ofexperimentally observed phenomena. One challenge of modeling such systems is realistic de-scription of the steric repulsion between nearby polymers. Previous models of the polymericexcluded volume interaction have the deficit of being non-analytic, being computationallyexpensive, or allowing polymers to erroneously cross each other. A recent solution to theseissues, implemented in the MEDYAN simulation platform, uses analytical expressions ob-tained from integrating an interaction kernel along the lengths of two polymer segments todescribe their repulsion. Here, we extend this model by re-deriving it for lower-dimensionalgeometrical configurations, deriving similar expressions using a steeper interaction kernel,comparing it to other commonly used potentials, and showing how to parameterize thesemodels. We also generalize this new integrated style of potential by introducing a segmentalLennard-Jones potential, which enables modelling both attractive and repulsive interactionsin semi-flexible polymer networks. These results can be further generalized to facilitate thedevelopment of effective interaction potentials for other finite elements in simulations ofsoft-matter systems. a r X i v : . [ c ond - m a t . s o f t ] F e b raphical Abstract: The segmental Lennard-Jones potential allows for tunable attractive andrepulsive interactions between finite cylindrical segments.2 Introduction
Excluded volume interactions between spa-tially extended macromolecules play an impor-tant role in a wide range of cellular phenomena.They help to produce mesoscopically orderedstructures, which enable the complex function-ality exhibited by cells. For example, it hasbeen shown that steric interactions alone caninduce alignment of the long biopolymers thatcomprise the cytoskeleton [1]. These interac-tions have also been implicated in transmittingnon-equilibrium fluctuations from one cellularsubsystem to another [2]. Excluded volume,or steric, interactions are in fact an importantphysical feature of many soft matter systems,which are often controlled by a complex inter-play of steric and entropic effects [3, 4].The computational modeling of suchsoft matter systems has become an essen-tial tool used in biology, chemistry, andphysics [5, 6]. In particular, software pack-ages for simulating the cytoskeleton, a com-plex soft active matter system comprising in-terlinked biopolymers and molecular motors,have helped provide theoretical understand-ing of various experimental phenomena [7–12].Reaching timescales of thousands of secondsand length scales of tens of micrometers, pack-ages such as AFiNeS, CytoSim, the model ofKim and coworkers, and MEDYAN allow ex-ploration of fascinating emergent cytoskeletalphenomena while striving to preserve realisticmicroscopic physics [13–16]. These models typ-ically employ effective, coarse-grained poten-tials based on ideas from polymer physics. Forexample, treating a semi-flexible polymer (forwhich the typical polymer length is comparableto the persistence length) as a one-dimensionalpiecewise linear chain, the mechanical strainenergy can be straightforwardly decomposedinto stretching and bending terms which can becomputed using harmonic functions of the lin-ear segments’ positions [17]. It is less straight-forward, however, to model the potential en-ergy mediating the excluded volume interac-tion between neighboring polymers. This istreated in different ways between CytoSim, themodel of Kim and coworkers, and MEDYAN,and it is not considered in AFiNeS. However,accurately modelling repulsion between poly- mers is essential for realistically simulating im-portant behaviors such as entanglement, repta-tion, liquid crystal ordering, and entropic de-pletion forces [4, 18–20].The primary physical origin of theexcluded volume interaction between typi-cal biopolymers such as actin is screenedCoulomb repulsion [21–23]. Actin filamentshave a relatively high linear charge density( ∼ . e/ ˚ A ), but biological ionic environmentshave a Debye-H¨uckel screening length ( ∼ ∼ . ∼ . − ∼ − r can be used as a smoothmimic of the hard-wall interaction. Specifi-cally, in MEDYAN the function 1 /r servesas the interaction kernel of the double integralover the lengths of the two segments (see Equa-tion 1 below). However, the result of the inte-gration is an opaque and complicated expres-sion, and it contains degeneracies when the twolinear segments are coplanar, leading to unde-fined behavior. To address these shortcomings,in this paper we first clarify the calculation ofthe excluded volume repulsion potential usedin MEDYAN. Then we illustrate how the prob-lem can be solved in the coplanar case and inother lower-dimensional geometries, and derivefurther expressions for the alternative steeperinteraction kernel 1 /r . We then character-ize the dependence of these interactions on theconfigurations of the two segments, discuss howto parameterize the potential, and compare itto the widely-used Gay-Berne form. We alsointroduce a new ‘segmental Lennard-Jones’ in-teraction which has both attractive and repul-sive components. Finally, we implement a nu-merical approximation method and discuss thegain in computational efficiency from using theanalytical expressions. Here we derive analytical expressions for theexcluded volume repulsion energy between twopolymer segments, using the 1 /r interactionkernel which corresponds to the implementa-tion in MEDYAN. We first give the derivationfor cylindrical segments in 3D space. Becausethe 3D expressions for the repulsion energy arenot defined when the cylinders are coplanar,we next describe the steps for re-deriving theseexpressions in 2D scenarios. Finally, we extendthe derivations to apply to an interaction ker-nel of 1 /r , representing an even steeper hard-wall mimic. Throughout this section certaincomplicated integrals must be solved, for whichwe use the computer algebra system (CAS)Mathematica [33, 34]. We provide Mathemat-ica notebook (.nb) files in the Supplementary Material which implement the calculations de-scribed below. Given the positions and orientations of twothin cylinders, we define the excluded volumerepulsion energy U as proportional to a doubleintegral of the function 1 /r ( s, t ) , where theintegrals run over the length of each cylinder: U = K vol (cid:90) (cid:90) dsdt r ( s, t ) . (1)Here s , t ∈ [0 ,
1] parameterize the distancealong the two cylinders A and B respectively, r ( s, t ) denotes the magnitude of the vector r ( s, t ) which separates points on the two cylin-der axes p s and p t , and K vol is the constant ofproportionality, having units of energy timeslength squared. The geometry of the problemis illustrated in Figure 1.Figure 1: Two cylinders, A and B , are in red,with the minus ends marked by blue dots. r isdefined by the points p s and p t on the cylin-ders. The position and orientation of thetwo cylinders A and B can be specified by fourvectors pointing to the positions of the fourends: c and c point to the minus and plusends, respectively, of cylinder A and c and c do likewise for cylinder B . Alternatively, we4an describe the two cylinders by the vectors V = c (2) A = c − c (3) B = c − c (4) C = c − c . (5)We can represent a point on cylinder A param-eterized by s , p s , as p s = V − s A , (6)and similarly for a point on cylinder B param-eterized by t , p t , we have p t = V + C − t B . (7)To solve the integral in Equation 1, first weneed to write r ( s, t ). We have r = p t − p s = C − t B + s A (8)and r = (cid:0) C · C + 2 s A · C + s A · A − t B · C − st A · B + t B · B (cid:1) / . (9)To simplify notation, we introduce the follow-ing variables: a = A · A b = B · B c = C · C d = A · B e = A · C f = B · C . Defining these intermediary variables is alsocomputationally efficient by avoiding repeat-edly calculating the same expressions. Withthis, the goal is do the following integral: U = K vol (cid:90) (cid:90) dsdt c + 2 es + as − f t − dst + bt ) . (10)This integral can be done with the help ofa CAS resulting in a lengthy expression pro-vided in the Appendix (Equation 27). It can be cleaned up somewhat by introducing the fol-lowing variables: AA = (cid:112) ac − e BB = (cid:112) bc − f CC = de − afDD = be − dfEE = (cid:112) a ( b + c − f ) − ( d − e ) F F = (cid:112) b ( a + c + 2 e ) − ( d + f ) GG = d − ab − CCHH = CC + GG + DDJ J = c ( GG + CC ) + eDD − f CCAT G = tan − (cid:18) a + eAA (cid:19) − tan − (cid:16) eAA (cid:17) AT G = tan − (cid:18) a + e − dEE (cid:19) − tan − (cid:18) e − dEE (cid:19) AT G = tan − (cid:18) fBB (cid:19) − tan − (cid:18) f − bBB (cid:19) AT G = tan − (cid:18) d + fF F (cid:19) − tan − (cid:18) d + f − bF F (cid:19) . These variable names are chosen to matchthose in the MEDYAN codebase. We pointout that the multiple letters comprising thesevariables do not indicate multiplication of twoor more variables. With this, the result is U = K vol J J (cid:18)
AT G CCAA + AT G GGEE + AT G DDBB + AT G HHF F (cid:19) . (11)Implementing excluded volume repul-sion in simulation usually also requires expres-sions for the derivatives of the energy withrespect to the cylinder endpoints c i , whichare used to determine the forces for timeintegrator-based approaches (e.g. CytoSim)or equivalently the gradients for minimization-based approaches (e.g. MEDYAN). Deriva-tives such as ∂U∂ c can be found using the chainrule: ∂U∂ c = ∂U∂a ∂a∂ c + ∂U∂b ∂b∂ c + . . . , (12)where the derivatives ∂U∂a , ∂U∂b , . . . , can be ob-tained from Equation 27 in the Appendix, and5he derivatives ∂a∂ c i can be found using the def-initions given above. For example, ∂a∂ c = ∂∂ c ( c · c − c · c + c · c )= c − c . (13) When the cylinders A and B are in the sameplane then the vectors A , B , and C are allcoplanar and the problem becomes effectively2D. Some implementations may also simply as-sume a 2D space. In this scenario, the scalartriple product ( A × B ) · C vanishes. It can beshown by straightforward algebraic rearrange-ment that the quantity J J appearing in thedenominator of the right hand side of Equa-tion 11 is given by
J J = − (( A × B ) · C ) . (14)As a result, the above expressions for the en-ergy U and derivatives ∂U∂ c i are not defined, anda special case must be considered.It is instructive to count the numberof free variables in the 3D and 2D case. In bothsettings, U is invariant with respect to a rigidrotation or translation of the system. In 3D,we originally have 12 variables (the 12 com-ponents of c , c , c , and c ), but translationinvariance implies that 3 degrees of freedom areextraneous and rotation invaraince implies that3 additional degrees of freedom are extraneous.This leaves 6 independent degrees of freedom,which appear in the expression for U as a , b , c , d , e , f . In 2D, we originally have 8 degreesof freedom, but translation invariance impliesthat 2 degrees of freedom are extraneous androtation invariance implies than additional 1degree of freedom is extraneous, leaving 5 de-grees of freedom. Indeed, the condition in 2Dthat J J = 0 implies an additional constraintamong the 6 variables. It can be shown that,in the 2D case but not in the 3D case, f = de + (cid:112) ( ab − d )( ac − e ) a , (15)and hence only 5 variables are free in 2D.One could through substitution writethe integrand 1 /r in terms of five free vari-ables in 2D, but this becomes an algebraically complicated expression that precludes exact in-tegration. Instead, we first rotate the con-figuration so the shared plane coincides withthe xy plane. Next, we write the integrandusing the 6 (redundant) vector components A x , A y , B x , B y , C x , and C y to find the inter-action energy in 2D. We have U = K vol (cid:90) (cid:90) dsdt C x + A x s − B x t ) + ( C y + A y s − B y t ) ) . (16)This integral has a complicated result which isprovided in the Appendix (Equation 28). Thedenominator of the result is proportional to A y B x − A x B y . If A and B are parallel (oranti-parallel) in addition to coplanar, then onecan show that A x B y = A y B x , and hence theexpression for U in the coplanar case is not de-fined. When A and B are (anti-)parallel,then B = ξ A for some ξ ∈ IR, ξ (cid:54) = 0. Ex-pressing the integrand using this new variable,we have U = K vol (cid:90) (cid:90) dsdt C x + A x ( s − ξt )) + ( C y + A y ( s − ξt )) ) . (17)The result of this integral is also provided inthe Appendix (Equation 29). The denomina-tor of that result is proportional to ( A y C x − A x C y ) . If, in addition to being parallel, A and B are colinear, then C is parallel to A and A x C y = A y C x , and this result is not defined.When A and B are colinear, one maywrite C = ζ A and express the integrand as U = K vol (cid:90) (cid:90) dsdt A x ( ζ + s − ξt )) + ( A y ( ζ + s − ξt )) ) . (18)The result of this integral, also provided in theAppendix (Equation 30), is simpler than in theprevious cases, depending on just 4 variables A x , A y , ξ, and ζ . Several ratios appear inthe result with denominators proportional to( A x + A y ) , ξ, ζ, ( ζ − ξ ) , ( ζ − ξ + 1) , and( ζ + 1) . Requiring that none of these are zero6mplies that A , B , and C are all non-zero, andthat A and B are nowhere coincident in whichcase the interaction would diverge.The treatment given above for han-dling the special case scenarios of cylinder con-figurations that lead to degeneracies in the en-ergy expressions is not exhaustive, and certaindegeneracies remain (such as one coming from A x C x + A y C y = 0 in the parallel case, Equation29). These degeneracies, unlike those comingfrom geometrical properties of the configura-tions, remain as a result of expressing the in-tegrand using redundant variables (which wasnecessary to do the integration analytically).It would be straightforward to consider eachdegeneracy in turn and, by following steps sim-ilar to those outlined above, derive backupexpressions for each scenario. However, werecommend instead implementing a numericalapproximation method to fall back on whenthese degeneracies are encountered in simula-tion. Such a numerical method is described be-low. We emphasize that in a given dimension-ality, the manifold of cylinder configurationsleading to degeneracies is of lower dimensionthan the ambient space and hence such config-urations will be exceedingly rare under typicalphysical dynamics. /r interaction kernel In Equation 1, the repulsion energy betweentwo cylinders was taken as a double integralover both cylinder lengths of the interactionkernel 1 /r . The 1 /r interaction is fairlysteep, mimicking a hard-wall boundary withan effective cylinder radius set by the choiceof prefactor K vol . However one may prefer aneven steeper potential than 1 /r , such as 1 /r ,so that the range of separation over which theinteraction starts to be felt is narrower. Thenew interaction energy is expressed as U = K vol (cid:90) (cid:90) dsdt r ( s, t ) . (19)Such a potential mimics even more closely atrue hard-wall interaction (see Figure 3 below),and with both interactions in hand it becomespossible by combining them to create bimodalenergy profiles, similar to a Lennard-Jones po-tential. It is straightforward to carry throughidentical steps for the 1 /r kernel as outlined above for the 1 /r kernel, with the same is-sues of degeneracies arising from special-casecylinder configurations. For brevity, and sinceno new concepts are involved, we skip the dis-cussion here of how those steps are carried outand also omit the resulting expressions fromthe Appendix. The expressions can be foundin the supplementary Mathematica notebookfiles. Here we analyze the ‘integrated kernel’ energyfunctions (Equations 1 and 19) and discuss no-table features arising from a set of test cases.For comparison, we also introduce two other‘endpoint-based’ interaction functions which,rather than integrating the kernels 1 /r or 1 /r over the lengths of the cylinders, simply includerepulsion between felt by the endpoints of thetwo cylinders: U = K vol (cid:18) m ( c , c , c ) + 1 m ( c , c , c ) +1 m ( c , c , c ) + 1 m ( c , c , c ) (cid:19) (20)and U = K vol (cid:18) m ( c , c , c ) + 1 m ( c , c , c ) +1 m ( c , c , c ) + 1 m ( c , c , c ) (cid:19) (21)where m ( c x , c a , c b ) represents the minimal dis-tance from the point c x to the line segmentconnecting c a and c b . Expressions of this typeare sometimes used to model the steric repul-sion of polymers, but we show below that theyhave the deficiency of a relatively flat energyprofile for cylinder separations much less thanthe cylinder length, which can allow cylindersto overlap each other under typical dynam-ics [35, 36].We consider two cylinders each oflength 200, where the units are fixed by set-ting K vol = 1 throughout. The cylinders are7ligned (i.e. A · C = 0), and the vector l joiningthe each cylinder’s midpoint is kept perpendic-ular to each cylinder as we rotate one cylinderabout this vector, producing different relativeconfigurations. This set-up is visualized in theinsets of Figure 2.A-C. For each configurationwe vary the distance l (the magnitude of l ), andstudy the effect on the various interaction en-ergies Equations 1, 19, 20, and 21. The resultsare displayed in Figure 2.Two key features are evident fromthis example. First, there are different asymp-totic behaviors in the small and large distanceregimes, with a crossover around distances onthe order of L . For l (cid:28) L , the integrated ex-pressions behave like a power law with an expo-nent equal to that of the kernel function plus 2,whereas for l (cid:29) L , the integrated expressionsbehave like a power law with an exponent equalto that of the kernel function: U = (cid:90) (cid:90) dsdt r ( s, t ) n ∼ (cid:40) l − n l (cid:28) Ll − n l (cid:29) L (22)where n = 4 and 6 in the examples shown.This behavior is expected, since in the far fieldall points in the cylinders repel each otherwith similar magnitudes, whereas in the nearfield the repulsion is dominated only by nearbypoints, changing the scaling by a factor of l .Second, we see that for the non-coplanar cylin-der configurations (Figures 2.A and 2.B), theendpoint-based interaction energies have flatenergy profiles for distances much less than thecylinder length. For the endpoint-based func-tions, U ∼ (cid:40) l (cid:28) Ll − n l (cid:29) L . (23) One can understand this as resulting from thefact that, when l (cid:28) L , the distance betweenthe endpoints of non-coplanar cylinders changemuch less than the distance between pointsthe middle of the cylinder as l is decreased,and these points in the middle contribute donot contribute to the energy penalty in theendpoint-based case. This qualitative differ-ence between the integrated and endpoint-based interactions is much less pronounced inthe coplanar case, when the distance betweenthe endpoints change at the same rate as allpoints when l is decreased, as shown in Figure2.C. We next introduce a ‘segmentalLennard-Jones’ interaction potential U = K vol,6 (cid:18)(cid:90) (cid:90) dsdt r ( s, t ) (cid:19) n − K vol,4 (cid:18)(cid:90) (cid:90) dsdt r ( s, t ) (cid:19) n . (24)This expression has qualitative similarity tothe familiar 6-12 Lennard-Jones potential be-tween two particles, and can be tuned bychoosing the four parameters K vol,6 , K vol,4 , n , and n to mimic interactions that haveboth attractive and repulsive parts. For in-stance, computational modeling of depletionforces, which tend to aggregate polymers to-gether, may make use of an effective attrac-tive component in the polymer-polymer in-teraction [36]. An example of a segmentalLennard-Jones potential is illustrated in Fig-ure 2.D for the choices K vol,6 = K vol,4 = 1 and n = n = 1 / l is varied fordifferent cylinder orientations and choices of interaction energy. Units are arbitrary, as lengthand energy scales are set by the K vol prefactor which is taken to be unity (amounting to settingthe vertical position of each curve on the log-log plots). (A) The blue line represents Equation1, the red line represents Equation 19, the green line represents Equation 20, and the purpleline represents Equation 21. The dotted vertical line indicates the length of the two cylinders.The cylinders are oriented perpendicularly to each other as visualized in the inset. The lightdotted lines proportional to 1 /l , 1 /l , and 1 /l show the scaling behaviors in different regimes.(B) The same plot is shown as in panel (A), except for a relative orientation of 30 ◦ betweenthe cylinders. (C) The same plot is shown as in panel (A), except for a relative orientationof 0 ◦ between the cylinders. (D) A plot is shown of a segmental Lennard-Jones interactionenergy profile from combining the 1 /r and 1 /r integrated kernel energies (Equation 24) . Thehorizontal dotted line separates the attractive and repulsive regions. K vol Actual biopolymers can differ significantly intheir diameters, requiring that K vol be tunedfor particular biopolymers. This choice can bemade so that at the effective diameter d ∗ thetypical interaction energy U t (i.e. the energyfor some typical configuration of segments) isequal to some energetic penalty for steric over- lap U m of the system: U t ( d ∗ ) = U m . (25)Here we indicate how to use Equation 25 to de-termine K vol for actin and microtubules whichhave been modeled as chains of 100 nm longlinear segments. The radius of an actin fila-ment is approximately 3 . nm , and for a mi-9rotubule it is 12 . ◦ with respectto each other. In Figure 3 we show the resultof using U t ( d ∗ ) = U m to determine K vol , tak-ing U m = 41 pN nm to be 10 times the thermalenergy (as k B T = 4 . K vol = 4 . × pN nm for the 1 /r kerneland K vol = 4 . × pN nm for the 1 /r ker-nel. For microtubule segments, this proceduregives K vol = 7 . × pN nm for the 1 /r kernel and K vol = 7 . × pN nm for the1 /r kernel. We note that this parameteriza-tion process also depends on the chosen lengthof the cylindrical segments.Figure 3: Plots of the energy as a function ofthe separation distance are shown for differentparameterizations and interaction kernels. Thered and blue curves are parameterized for actinusing the 1 /r and 1 /r kernel respectively,while the green and purple curves are likewiseparamaterized for microtubules. The dashedvertical lines are drawn at d ∗ for actin and mi-crotubules, and the intersection points indicatethe enforcement of the condition U t ( d ∗ ) = U m where U m = 41 pN nm (horizontal dahsedline). K vol to typical con-figurations The parameterization method described abovehas one seemingly major ambiguity, which is how to determine the typical configuration ofthe segments at which to evaluate U t ( d ∗ ). Fix-ing the position of one segment and both seg-ments’ lengths, 5 variables remain to specifythe other segment: the offset vector C and thespherical coordinates θ (inclination) and φ (az-imuthal) of the unit vector ˆ B . The dependenceon the separation of aligned cylinders has beendiscussed above. We next explore the orienta-tional coordinates θ and φ , setting c = ( − L/ , − l, c = ( L/ , − l, c = (0 , , c = ( L sin θ cos φ, L sin θ sin φ, L cos θ ) . The geometry of this set up is illustrated in theinset of Figure 4.Figure 4: A surface plot of U as θ and φ arevaried is shown for the conditions described inthe main text. The colors indicate the energy,ranging from 0 pN nm (purple) to 50 pN nm(red), the gray region corresponds to configu-rations in which the segments nearly overlap,and the inset illustrates the definition of thespherical coordinates θ and φ .In Figure 4, we show the energy U for L = 100 nm, l = 10 nm, and K vol = 4 . × pN nm as a function of θ and φ . Evidently, U lies within a fairly small range over mostof the domain of θ and φ (with the exceptionbeing where the cylinder nearly overlap). It10an similarly be shown that the degree of free-dom corresponding to sliding one cylinder ina direction parallel to the other cylinder onlyaffects the interaction energy to within an or-der of magnitude. Therefore, we may concludethat the parameterization is not very sensitiveto how the typical configuration is chosen, andthat reliable order of magnitude estimates of K vol can be obtained for a given d ∗ and U m . Next, we compare the new integrated kernelexpression for the interaction energy betweencylindrical objects, Equation 1, to the widelyused Gay-Berne potential which describes theinteraction between anisotropic ellipsoidal ob-jects. The original Gay-Berne potential wasdesigned to be similar to a Lennard-Jonespotential, having both attractive and repul-sive contributions, but for comparison here wemodify the original potential to be only repul-sive and with an exponent of −
4; we give theformula for the Gay-Berne potential used herein the Appendix (Equation 31). We comparethe interaction energy profiles as a function ofdistance for two offset, rotated cylindrical seg-ments with variable aspect ratios. Holding thediameters d ∗ = 1 fixed, we change the cylin-der lengths L = κd ∗ , where κ is the geometricaspect ratio, and show that for large κ the Gay-Berne potential deviates strongly from the de-sired power-law repulsion. We use the follow-ing test case configuration, illustrated in theinset of Figure 5: for each choice of L , the hor-izontal offset (along their lengths) of two par-allel, initially aligned cylinders is chosen suchthat half of their lengths overlap, and one cylin-der is then rotated 45 ◦ around around the linejoining its midpoint and the other cylinder’soverlapping endpoint. The length of this line isthen varied to construct the interaction energyprofile for this test configuration. The energyscale is fixed by setting each energy to 1 at adistance of d ∗ . We display the results in Figure5. We observe that for large values of κ ,the Gay-Berne interaction profile deviates sig-nificantly from the expected power-law behav-ior, exhibiting weakened repulsion for l (cid:38) d ∗ and enhanced repulsion for l (cid:46) d ∗ compared to the integrated kernel interaction. On the otherhand, for κ = 1 the Gay-Berne profile and theintegrated kernel profile nearly coincide.Figure 5: Interaction energy profiles for theintegrated 1 /r kernel potential (top panel,Equation 1) and Gay-Berne potential (bottompanel, Equation 31), and the are shown. Inboth panels, the geometric aspect ratio κ isvaried in increments of 5 from 1 to 31 (i.e. L is increased with d ∗ fixed), as the colors arevaried from light blue to dark purple. The in-set shows the set-up of the two cylinders for κ = 21. The need for a new potential to de-scribe polymer repulsion can be understood asarising from the fact that, when modeling con-secutive cylinders in a polymer as ellipsoids,the potential energy is not uniform along thepolymer’s length. One can imagine a chain ofsausage links to represent this scenario. In the11ntegrated kernel interaction however, the en-ergy is uniform and therefore does not dependon how the polymer is discretized into cylin-ders. This is an important physical feature topreserve in computational modeling.Figure 6: The ratio of the numerically obtainedenergy U to the analytical U ∗ and the ratio ofthe numerical evaluation time t to that of theanalytical result t ∗ are shown as the number ofsampling points N sample is varied. The greenshaded area indicates where the agreement be-tween U and U ∗ is acceptable (i.e. U/U ∗ ≈ Finally, we illustrate the gain in computa-tional efficiency from having an analytical re-sult (Equation 11) for the integral in Equation1 rather than a numerical approximation, asis sometimes used in LAMMPS [32]. We im-plemented both the analytical result and a nu-merical scheme sampling the double integral at N sample points along each cylinder in compiledC code. The numerical scheme approximates U as U ≈ K vol N N sample (cid:88) i,j =1 r (cid:16) iN sample , jN sample (cid:17) , (26)where r ( s, t ) is given in Equation 9. For a sin-gle test case of aligned cylinders rotated by 45 ◦ relative to each other, we compared the energyand evaluation time for the numerical schemeto the analytical counterpart as N sample wasvaried from 5 to 30. The result is displayed inFigure 6. Once the number of sampling pointsis large enough that the numerical approxima-tion is acceptable ( N sample ∼ c i can be brought inside thesum in Equation 26. Our goal has been to clarify the derivation ofthe novel excluded volume repulsion potentialimplemented in MEDYAN and to extend thederivation to other scenarios of interest. Thisoverall approach to modeling repulsion interac-tions based on integrating an interaction ker-nel may be extended to other geometrical el-ements of finite size, such as 2D faces or 3Dvolumes. Despite the complexity of the re-sulting expressions for the energy and forces,they have the significant benefit of being an-alytical and avoiding endpoint-based interac-tions, which have flat energy profiles that canallow the repelling objects to erroneously over-lap each other. On the other hand, these ex-pressions have the issue of being undefined forcertain lower-dimensional rare configurationswhich impedes their usability in simulation.However, we have shown how this issue can behandled by re-deriving expressions using a re-duced number of variables. In addition, we de-scribed how other types of interactions can bedesigned, such as steeper repulsion and a seg-mental Lennard-Jones interaction, while stillaccounting for the finite dimensions of the in-teracting objects. This potential could be use-ful to model certain aggregating polymer sys-tems such as toroidal DNA [37,38]. The mathe-matical elaborations presented here should en-able other investigators to effectively use thesenew potentials in their computational studiesof soft matter systems.12 cknowledgements
We thank Wonyeong Jung and Tae YoonKim for helpful discussions. This work wassupported by the grants 1632976 and CHE-1800418 from the National Science Foundation.13 ppendix
Here we provide the full expressions of the interaction energy using the 1 /r interaction kernelin various types of cylinder configurations. The meaning of the variables is provided in the maintext. Cylinders in 3D U = K vol a ( f − bc ) + e ( be − df ) + cd ) (cid:32) (cid:0) ab + be − d − df (cid:1) tan − (cid:18) − d − f √ ab + b ( c +2 e ) − ( d + f ) (cid:19)(cid:112) ab + b ( c + 2 e ) − ( d + f ) − (cid:0) ab + be − d − df (cid:1) tan − (cid:18) b − d − f √ ab + b ( c +2 e ) − ( d + f ) (cid:19)(cid:112) ab + b ( c + 2 e ) − ( d + f ) + ( a ( b − f ) + d ( e − d )) tan − (cid:18) e − d √ a ( b + c − f ) − ( d − e ) (cid:19)(cid:112) a ( b + c − f ) − ( d − e ) + ( a ( f − b ) + d ( d − e )) tan − (cid:18) a − d + e √ a ( b + c − f ) − ( d − e ) (cid:19)(cid:112) a ( b + c − f ) − ( d − e ) + tan − (cid:16) e √ ac − e (cid:17) ( af − de ) √ ac − e + tan − (cid:16) a + e √ ac − e (cid:17) ( de − af ) √ ac − e + tan − (cid:18) b − f √ bc − f (cid:19) ( be − df ) (cid:112) bc − f + tan − (cid:18) f √ bc − f (cid:19) ( be − df ) (cid:112) bc − f (cid:33) (27) Coplanar cylinders U = K vol A y B x − A x B y ) (cid:32) ( A y B x − A x B y ) ( A y ( C x − B x ) + A x ( B y − C y )) ( A y C x − A x C y ) ( B y C x − B x C y )+ ( A y B x − A x B y ) ( A y ( C x − B x ) + A x ( B y − C y )) ( A y C x − A x C y ) ( A y B x + C y B x − A x B y − B y C x )+ tan − (cid:16) A x C x + A y C y A y C x − A x C y (cid:17) (cid:0) A x + A y (cid:1) ( A y C x − A x C y ) − tan − (cid:16) A x + C x A x + A y ( A y + C y ) A y C x − A x C y (cid:17) (cid:0) A x + A y (cid:1) ( A y C x − A x C y ) − tan − (cid:16) B x C x + B y C y B y C x − B x C y (cid:17) (cid:0) B x + B y (cid:1) ( B y C x − B x C y ) + tan − (cid:16) − B x + C x B x + B y ( C y − B y ) B y C x − B x C y (cid:17) (cid:0) B x + B y (cid:1) ( B y C x − B x C y ) + tan − (cid:16) A x B x + C x B x + B y ( A y + C y ) − A y B x − C y B x + A x B y + B y C x (cid:17) (cid:0) B x + B y (cid:1) ( A y B x + C y B x − A x B y − B y C x ) − tan − (cid:16) − B x + A x B x + C x B x + B y ( A y − B y + C y ) − A y B x − C y B x + A x B y + B y C x (cid:17) (cid:0) B x + B y (cid:1) ( A y B x + C y B x − A x B y − B y C x ) + tan − (cid:16) A x ( C x − B x )+ A y ( C y − B y ) A y ( B x − C x )+ A x ( C y − B y ) (cid:17) (cid:0) A x + A y (cid:1) ( A y ( B x − C x ) + A x ( C y − B y )) − tan − (cid:16) A x +( C x − B x ) A x + A y ( A y − B y + C y ) A y ( B x − C x )+ A x ( C y − B y ) (cid:17) (cid:0) A x + A y (cid:1) ( A y ( B x − C x ) + A x ( C y − B y )) (cid:33) (28)14 arallel cylinders U = K vol ξ ( A y C x − A x C y ) (cid:32) A x (cid:32) tan − (cid:18) A x + C x A x + A y ( A y + C y ) A y C x − A x C y (cid:19) + ξ tan − (cid:18) − ξA x + C x A x + A y ( C y − ξA y ) A x C y − A y C x (cid:19) + ξ tan − (cid:18) (1 − ξ ) A x + C x A x + A y ( − ξA y + A y + C y ) A y C x − A x C y (cid:19) + tan − (cid:18) (1 − ξ ) A x + C x A x + A y ( − ξA y + A y + C y ) A x C y − A y C x (cid:19) (cid:19)(cid:33) + C x A x (cid:32) tan − (cid:18) A y C x − A x C y A x C x + A y C y (cid:19) − tan − (cid:18) A y C x − A x C y A x + C x A x + A y ( A y + C y ) (cid:19) + tan − (cid:18) − ξA x + C x A x + A y ( C y − ξA y ) A y C x − A x C y (cid:19) − tan − (cid:18) (1 − ξ ) A x + C x A x + A y ( − ξA y + A y + C y ) A y C x − A x C y (cid:19) (cid:19)(cid:33) + A y (cid:32) A y (cid:32) tan − (cid:18) A x + C x A x + A y ( A y + C y ) A y C x − A x C y (cid:19) + ξ tan − (cid:18) − ξA x + C x A x + A y ( C y − ξA y ) A x C y − A y C x (cid:19) + ξ tan − (cid:18) (1 − ξ ) A x + C x A x + A y ( − ξA y + A y + C y ) A y C x − A x C y (cid:19) + tan − (cid:18) (1 − ξ ) A x + C x A x + A y ( − ξA y + A y + C y ) A x C y − A y C x (cid:19) (cid:19)(cid:33) + C y (cid:18) tan − (cid:18) A y C x − A x C y A x C x + A y C y (cid:19) − tan − (cid:18) A y C x − A x C y A x + C x A x + A y ( A y + C y ) (cid:19) + tan − (cid:18) − ξA x + C x A x + A y ( C y − ξA y ) A y C x − A x C y (cid:19) − tan − (cid:18) (1 − ξ ) A x + C x A x + A y ( − ξA y + A y + C y ) A y C x − A x C y (cid:19) (cid:33)(cid:33)(cid:33) (29) Colinear cylinders U = K vol ξ (cid:0) A x + A y (cid:1) (cid:32) − ζ + 1( ζ − ξ ) − ζ − ξ + 1) + 1( ζ + 1) (cid:33) (30) The Gay-Berne potential
The Gay-Berne potential is designed to generalize the familiar Lennard-Jones interaction togeometrically anisotropic ellipsoidal particles. Further generalizations to lower symmetry inter-actions have also been constructed, but we assume here a pair of identical radially symmetricellipsoids repelling with a 1 /r potential. The formulas given here are adapted from Ref 39.The interaction energy is written as U ( ˆu i , ˆu j , r ij ) = 4 (cid:15) (cid:15) ( ˆu i , ˆu j , ˆr ij ) (cid:18) σ s r ij − σ ( ˆu i , ˆu j , ˆr ij ) + σ s (cid:19) . (31)Here r ij points from the center of ellipsoid i to the center of ellipsoid j , u i points along themajor axis of ellipsoid i and likewise for u j , the caret hats indicate unit vectors, and r ij is themagnitude of r ij . σ s represents the length of the minor ellipsoid axis (the ‘diameter’), and σ e represents the length of the major axis (the ‘length’). The prefactor (cid:15) sets the energy scale.The shape function σ ( ˆu i , ˆu j , ˆr ij ) is σ ( ˆu i , ˆu j , ˆr ij ) = σ s (cid:18) − χ (cid:18) ( ˆr ij · ˆu i + ˆr ij · ˆu j ) χ ˆu i · ˆu j + ( ˆr ij · ˆu i − ˆr ij · ˆu j ) − χ ˆu i · ˆu j (cid:19)(cid:19) − / , (32)15here χ = κ − κ +1 and κ = σ e σ s . The interaction function (cid:15) ( ˆu i , ˆu j , ˆr ij ) is (cid:15) ( ˆu i , ˆu j , ˆr ij ) = (cid:0) − χ ( ˆu i · ˆu j ) (cid:1) − ν (cid:18) − χ (cid:48) (cid:18) ( ˆr ij · ˆu i + ˆr ij · ˆu j ) χ (cid:48) ˆu i · ˆu j + ( ˆr ij · ˆu i − ˆr ij · ˆu j ) − χ (cid:48) ˆu i · ˆu j (cid:19)(cid:19) µ , (33)where χ (cid:48) = k (cid:48) /µ − k (cid:48) /µ +1 , k (cid:48) = (cid:15) s (cid:15) e , and (cid:15) s and (cid:15) e represent, respectively, the depth of the potential wellfor the side-to-side and end-to-end configurations of the two ellipsoids. The free parameters ofthis energy are (cid:15) , (cid:15) s , (cid:15) e , σ s , σ e , and the fitting exponents µ and ν . For the comparisons donein Figure 5, we take σ s = 1, (cid:15) e = (cid:15) s = 1, µ = 2 and ν = 1 (following Ref. 39), and (cid:15) is chosenso that U = 1 when the separation is d ∗ , as described in the main text. σ e = κσ s is varied totest the effect of geometrical anisotropy. References [1] Magdalena Preciado L´opez, Florian Huber, Ilya Grigoriev, Michel O Steinmetz, AnnaAkhmanova, Marileen Dogterom, and Gijsje H Koenderink. In vitro reconstitution ofdynamic microtubules interacting with actin filament networks. In
Methods in Enzymology ,volume 540, pages 301–320. Elsevier, 2014.[2] Clifford P Brangwynne, Gijsje H Koenderink, Frederick C MacKintosh, and David A Weitz.Nonequilibrium microtubule fluctuations in a model cytoskeleton.
Physical Review Letters ,100(11):118104, 2008.[3] Masao Doi.
Soft Matter Physics . Oxford University Press, 2013.[4] Michael Rubinstein, Ralph H Colby, et al.
Polymer Physics , volume 23. Oxford universitypress New York, 2003.[5] Ulf D Schiller, Timm Kr¨uger, and Oliver Henrich. Mesoscopic modelling and simulation ofsoft matter.
Soft Matter , 14(1):9–26, 2018.[6] Thomas E Gartner III and Arthi Jayaraman. Modeling and simulations of polymers: Aroadmap.
Macromolecules , 52(3):755–786, 2019.[7] Jonathon Howard et al. Mechanics of motor proteins and the cytoskeleton. 2001.[8] Xiaona Li, Qin Ni, Xiuxiu He, Jun Kong, Soon-Mi Lim, Garegin A Papoian, Jerome P Trze-ciakowski, Andreea Trache, and Yi Jiang. Tensile force-induced cytoskeletal remodeling:Mechanics before chemistry.
PLoS Computational Biology , 16(6):e1007693, 2020.[9] Aravind Chandrasekaran, Arpita Upadhyaya, and Garegin A Papoian. Remarkable struc-tural transformations of actin bundles are driven by their initial polarity, motor activity,crosslinking, and filament treadmilling.
PLoS Computational Biology , 15(7):e1007156, 2019.[10] Carlos Floyd, Garegin A Papoian, and Christopher Jarzynski. Quantifying dissipation inactomyosin networks.
Interface Focus , 9(3):20180078, 2019.[11] Simon L Freedman, Cristian Suarez, Jonathan D Winkelman, David R Kovar, Gregory AVoth, Aaron R Dinner, and Glen M Hocky. Mechanical and kinetic factors drive sort-ing of f-actin cross-linkers on bundles.
Proceedings of the National Academy of Sciences ,116(33):16192–16197, 2019.[12] Julio M Belmonte, Maria Leptin, and Fran¸cois N´ed´elec. A theory that predicts behaviorsof disordered cytoskeletal networks.
Molecular Systems Biology , 13(9):941, 2017.1613] Simon L Freedman, Shiladitya Banerjee, Glen M Hocky, and Aaron R Dinner. A versa-tile framework for simulating the dynamic mechanical structure of cytoskeletal networks.
Biophysical Journal , 113(2):448–460, 2017.[14] Francois Nedelec and Dietrich Foethke. Collective langevin dynamics of flexible cytoskeletalfibers.
New Journal of Physics , 9(11):427, 2007.[15] Taeyoon Kim, Wonmuk Hwang, Hyungsuk Lee, and Roger D Kamm. Computational anal-ysis of viscoelastic properties of crosslinked actin networks.
PLoS Computational Biology ,5(7):e1000439, 2009.[16] Konstantin Popov, James Komianos, and Garegin A Papoian. Medyan: Mechanochemicalsimulations of contraction and polarity alignment in actomyosin networks.
PLoS Compu-tational Biology , 12(4):e1004877, 2016.[17] Chase P Broedersz and Fred C MacKintosh. Modeling semiflexible polymer networks.
Reviews of Modern Physics , 86(3):995, 2014.[18] Pierre-Giles de Gennes. Reptation of a polymer chain in the presence of fixed obstacles.
The Journal of Chemical Physics , 55(2):572–579, 1971.[19] B G¨otzelmann, Robert Evans, and Siegfried Dietrich. Depletion forces in fluids.
PhysicalReview E , 57(6):6785, 1998.[20] Paul M Chaikin, Tom C Lubensky, and Thomas A Witten.
Principles of condensed matterphysics , volume 10. Cambridge university press Cambridge, 1995.[21] Rob Phillips, Jane Kondev, Julie Theriot, and Hernan Garcia.
Physical biology of the cell .Garland Science, 2012.[22] Michael Rubinstein and Garegin A Papoian. Polyelectrolytes in biology and soft matter.
Soft Matter , 8(36):9265–9267, 2012.[23] Rafael Tadmor, Ernesto Hern´andez-Zapata, Nianhuan Chen, Philip Pincus, and Jacob NIsraelachvili. Debye length and double-layer forces in polyelectrolyte solutions.
Macro-molecules , 35(6):2380–2388, 2002.[24] Thomas E Angelini, Ramin Golestanian, Robert H Coridan, John C Butler, Alexandre Be-raud, Michael Krisch, Harald Sinn, Kenneth S Schweizer, and Gerard CL Wong. Counteri-ons between charged polymers exhibit liquid-like organization and dynamics.
Proceedingsof the National Academy of Sciences , 103(21):7962–7967, 2006.[25] DA Smith and DG Stephenson. An electrostatic model with weak actin-myosin attach-ment resolves problems with the lattice stability of skeletal muscle.
Biophysical Journal ,100(11):2688–2697, 2011.[26] Fr´ed´eric Eghiaian, Annafrancesca Rigato, and Simon Scheuring. Structural, mechani-cal, and dynamical variability of the actin cortex in living cells.
Biophysical Journal ,108(6):1330–1340, 2015.[27] Lori K Sanders, Camilo Gu´aqueta, Thomas E Angelini, Jae-Wook Lee, Scott C Slimmer,Erik Luijten, and Gerard CL Wong. Structure and stability of self-assembled actin-lysozymecomplexes in salty water.
Physical Review Letters , 95(10):108302, 2005.[28] Paul A Janmey, Soren Hvidt, J K¨as, Dietmar Lerche, Anthony Maggs, Erich Sackmann,Manfred Schliwa, and Thomas P Stossel. The mechanical properties of actin gels. elasticmodulus and filament motions.
Journal of Biological Chemistry , 269(51):32503–32513,1994. 1729] JG Gay and BJ Berne. Modification of the overlap potential to mimic a linear site–sitepotential.
The Journal of Chemical Physics , 74(6):3316–3319, 1981.[30] Roberto Berardi, Carlo Fava, and Claudio Zannoni. A gay–berne potential for dissimilarbiaxial particles.
Chemical Physics Letters , 297(1-2):8–14, 1998.[31] Timothy W Sirk, Yelena R Slizoberg, John K Brennan, Martin Lisal, and Jan W Andzelm.An enhanced entangled polymer model for dissipative particle dynamics.
The Journal ofChemical Physics , 136(13):134903, 2012.[32] Steve Plimpton. Fast parallel algorithms for short-range molecular dynamics.
Journal ofComputational Physics , 117(1):1–19, 1995.[33] Wolfram Research, Inc. Mathematica, Version 12.2. Champaign, IL, 2020.[34] Albert Rich, Patrick Scheibe, and Nasser M Abbasi. Rule-based integration: An extensivesystem of symbolic integration rules.
Journal of Open Source Software , 3(32):1073, 2018.[35] Jamie Rickman, Fran¸cois N´ed´elec, and Thomas Surrey. Effects of spatial dimensionality andsteric interactions on microtubule-motor self-organization.
Physical Biology , 16(4):046004,2019.[36] Ga¨elle Letort, Antonio Z Politi, Hajer Ennomani, Manuel Th´ery, Francois Nedelec, andLaurent Blanchoin. Geometrical and mechanical properties control actin filament organi-zation.
PLoS Computational Biology , 11(5):e1004245, 2015.[37] Victor A Bloomfield. Dna condensation.
Current Opinion in Structural Biology , 6(3):334–341, 1996.[38] Am´elie Leforestier and Fran¸coise Livolant. Structure of toroidal dna collapsed inside thephage capsid.
Proceedings of the National Academy of Sciences , 106(23):9157–9162, 2009.[39] Wenduo Chen, Youliang Zhu, Fengchao Cui, Lunyang Liu, Zhaoyan Sun, Jizhong Chen,and Yunqi Li. Gpu-accelerated molecular dynamics simulation to study liquid crystal phasetransition using coarse-grained gay-berne anisotropic potential.