A continuum theoretical model and finite elements simulation of bacterial flagellar filament phase transition

A continuum theoretical model and finite elements simulation of bacterial flagellar filament phase transition

Accepted Manuscript Review A Continuum Theoretical Model and Finite Elements Simulation of Bacterial Flagellar Filament Phase Transition Xiaoling Wang...

2MB Sizes 0 Downloads 41 Views

Accepted Manuscript Review A Continuum Theoretical Model and Finite Elements Simulation of Bacterial Flagellar Filament Phase Transition Xiaoling Wang, Shuo Meng, Jingshi Han PII: DOI: Reference:

S0021-9290(17)30468-2 http://dx.doi.org/10.1016/j.jbiomech.2017.09.012 BM 8371

To appear in:

Journal of Biomechanics

Received Date: Revised Date: Accepted Date:

16 February 2017 6 September 2017 6 September 2017

Please cite this article as: X. Wang, S. Meng, J. Han, A Continuum Theoretical Model and Finite Elements Simulation of Bacterial Flagellar Filament Phase Transition, Journal of Biomechanics (2017), doi: http://dx.doi.org/10.1016/ j.jbiomech.2017.09.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Continuum Theoretical Model and Finite Elements Simulation of Bacterial Flagellar Filament Phase Transition Xiaoling Wang1, 2 #,*, Shuo Meng1,*, Jingshi Han1 1

School of Mechanical Engineering, University of Science and Technology Beijing, China, Beijing

100083, China 2

School of Engineering and Applied Sciences, Harvard University, Cambridge MA 02138, USA

*Equal contributors #

Corresponding author: [email protected]

Abstract The Bacterial flagellar filament can undergo a polymorphic phase transition in response to both mechanical and chemical variations in vitro and in vivo environments. Under mechanical stimuli, such as viscous flow or forces induced by motor rotation, the filament changes its phase from left-handed normal (N) to right-handed semi-coiled (SC) via phase nucleation and growth. Our detailed mechanical analysis of existing experiments shows that both torque and bending moment contribute to the filament phase transition. In this paper, we establish a non-convex and non-local continuum model based on the Ginzburg-Landau theory to describe main characteristics of the filament phase transition such as new-phase nucleation, growth, propagation and the merging of neighboring interfaces. The finite element method (FEM) is adopted to simulate the phase transition under a displacement-controlled loading condition (rotation angle and bending deflection). We show that new-phase nucleation corresponds to the maximum torque and bending moment at the stuck end of the filament. The hysteresis loop in the loading and unloading curves indicates energy dissipation. When 1

the new phase grows and propagates, torque and bending moment remain static. We also find that there is a drop in load when the two interfaces merge, indicating a concomitant reduction in the interfacial energy. Finally, the interface thickness is governed by the coefficients of the gradient of order parameters in the non-local interface energy. Our continuum theory and the finite element method provide a method to study the mechanical behavior of such biomaterials. Keywords: Bacterial flagellar filament, Phase transition, Ginzburg-Landau theory, Finite element method 1. Introduction The flagellar movement of live bacteria was first discovered by Ehrenberg (Lazarus and Jahn, 1998) in 1838. Many types of bacteria, such as Escherichia coli, Bacillus subtilis or Salmonella, are externally flagellated, where each flagellum consists of a long helical filament that is attached to a rotary motor embedded in the rigid cell wall via a flexible hook. These bacteria propel themselves forward by rotating a bundle of helical filaments as shown in Fig. 1(a). The flagellar filament is a large homogeneous assembly of a single protein called flagellin (Wang, 2006; Wang and Sun, 2010b). Macroscopically the filament can be viewed as a hollow tube with an outer diameter of 20 nm and a length of 15-20 μm (Turner et al., 2000). A filament has twelve forms with different radii and pitch lengths, with the difference in their microstructures characterized by protein-subunit misfit (e.g. 0.8×10-10 m). When the motor of each

flagella is in counter clockwise (CCW) rotation, all flagella form a rotating bundle

and propel the cell in one direction; when even one motor’s rotation reverses, the flagella fall apart, the cell tumbles, and the filament undergoes polymorphic transition. Upon reversal of

2

motor rotation , all flagella re-bundle and the cell swims in another direction, (Asakura, 1969; Asakura and Iino, 1972; Calladine, 1976; Hotani, 1976; Kamiya and Asakura, 1976b; Kamiya and Asakura, 1976a; Calldine, 1978; Hotani, 1982b; Darnton and Berg, 2007; Darnton et al., 2007) as shown in Fig. 1(a). The phase transition of the filament can be induced by various stimuli such as a change in temperature or pH, or by external forces (Darnton et al., 2007). Hotani’s experiment (Hotani, 1982b) showed cyclic chirality transformations of flagella tethered at one end to a glass surface and subjected to an external fluid flow, as shown in Fig. 1(b) and 1(c). We obtained the main component forces at the cross along the filament using the experimental parameters in Ref (Hotani, 1982b), as shown in Fig. 2(c-f). Tension force F in Fig. 2(f) causes linear elastic deformation along the filament; and bending moment M2 in Fig. 2(e) straightens the filament. Torque T in Fig. 2 (c) and bending moment M1 in Fig. 2(d) contribute to the phase nucleation. The torque T and the bending moment M1 of the filament reach their peak values at the stuck end, causing a new form to nucleate at that end, as shown in each transition cycle. This suggests that the new phase nucleation at the stuck end corresponds to the accumulations of both torque T and bending M1 (Wang, 2006). Several theoretical models have been proposed to describe filament phase transition. Goldstein et al. (Goldstein et al., 2000) and Coombs et al. (Coombs et al., 2002) extended Kirchhoff’s classic theory of an elastic rod by introducing a double-well potential for the spontaneous torsion to describe the chirality transition of a filament observed in Hotani’s experiment qualitatively. Power et al. (Srigiriraju and Powers, 2005; Srigiriraju and Powers, 2006) proposed another continuum theory for the polymorphism of the filament to account

3

for the alternate stable conformations. Wada et al. (Wada and Netz, 2008) introduced a bi-stable helical filament model that accounts for polymorphic states by a discrete Ising-like spin variable along the arc-length, and implemented the model to simulate the phase transition by hybrid Brownian dynamics and Monte-Carlo methods. Wang et al. (Wang et al., 2011) introduced a theoretical model of the non-convex non-local viscoelasticity and implemented it in a one dimensional (1D) non-local finite element method (FEM) code to simulate the filament’s phase transition. Here we extend the latter work to two dimensions (2D) by considering two order parameters in free energy. We implement the 2D free energy theoretical framework in a 2D non-local FEM code through proportional loading of rotation angle and bending deflection (in the ratio of 0.57) to simulate the filament’s phase transition between a left-handed normal (N) phase and a right-handed semi-coiled (SC) phase. The theoretical model and the FEM Methodology are outlined in the next two sections. The paper concludes with a discussion of the FEM simulation results. 2. Theoretical Model The Lagrange formulation for non-local and non-convex free energy (Reid and Gooding, 1994; Curnoe and Jacobs, 2001; Sun and He, 2008; He and Sun, 2009b) is

U 0  qi

(3)

where U  (W  G)dx , x is the spatial coordinate;  qi represent the variations of the freedom; W is the nonconvex elastic energy density that describes material instability and phase transition (Cahn and Hilliard, 1958; Falk, 1983; Barsch and Krumhansl, 1984), G is the non-local gradient energy density to account for the interfacial energy (Cahn and Hilliard, 4

1958; Mindlin and Eshel, 1968; Falk, 1983; Barsch and Krumhansl, 1984). To use the framework to describe the flagellar phase transition, suitable energy formulations are needed—non-convex elastic energy and elastic gradient energy. We use the generalized Kirchhoff theory to describe filament deformation ; and the Landau (non-convex) free energy to describe the elastic energy of the filament. The analysis above suggests that both torque and bending moment contribute to the phase transition under viscous fluid flow along the bacterial filament; therefore, we choose macroscopic curvature and twist as the order parameters in the filament free energy construction.

2.1 The Kirchhoff Model

A filament is defined as a three dimensional body whose cross-section is much smaller than its length. In addition, if the curvature of the filament is assumed to be small relative to its length, it is possible to show that its dynamics is governed by the celebrated Kirchhoff equations. We choose the orthonormal triad

ei 

along the filament as follows. With e3 along the

unit tangent vector, e1 and e2 will complete the local orthonormal system, but will have a physical meaning according to the microscopic structure of the filament. We can make a painted stripe on the side of the filament, use it to indicate the center of the R-type proto-filaments, with e2 always pointing to this stripe during the transition process, and e1 completing a right-handed triad, and the unit tangent vectors e3, e1, e2, e3, together forming the orthonormal system. Then if a filament is transformed from the left handed normal form to the right handed semi-coiled form, e2 will always point to the center of the R-type

5

proto-filaments, and e1 will always be perpendicular to e2 in the plane that is orthonormal to e3. This is shown in Fig. 3. Now we define the vector Φ(s) based on the angles of rotation of the system (e1, e2, e3) around some fixed frame of reference; at the point r(s), s is the arc length of the un-deformed axis. Any smooth deformation of the filament is now described by

  s    1  s  ,  2  s  ,  3  s   

d ds

(4)

Ω measures the rate of rotation of the e frame:

ei    ei s

(5)

where, Ω3(s) is the rate of rotation of the stripe around the filament centerline, Ω 2(s) measures the curvature in the plane containing the stripe, and Ω1(s) is the curvature in the direction perpendicular to the stripe. We will henceforth denote Ω3 as the twist rate. This description of the filament is sufficient for any differentiable configuration r with differentiable twist information, and includes information about the physical twist state. The Frenet-Serret curvature of filament is described in Fig.4. Henceforth, we shall use κ to represent the curvature Ω1(s), and γ to represent the twist rate Ω3(s).

2.2 2D non-convex free energy

A Landau (non-convex) free energy is used to describe the elastic energy of the filament. Since, from the analysis of the experimental work of Hotani, M2 does not contribute the phase transition of the filament, we assume that the curvature Ω2 deforms elastically.  and  are respectively curvature due to M1 and twist induced by T, with both contributing to the filament phase transition. Therefore, the free energy of the filament for small deformations can be simplified and expressed as a Taylor expansion in Ωi, where the indices are 1 and 3.

6

W1  

L

0

1 1 1 1 ( A222 2  Ai i   Aij i  j   Aijk i  j k   Aijkl i jk l  )ds 2 2 i, j 6 i , j ,k 24 i , j ,k ,l

(6)

Here, we use κ represents curvature Ω1(s), γ represents twist Ω3(s). We expand W1 to the fourth order as follows,

W1  a2 2  s  t  b 2  c 2  d  e 3  f  3  g 2  h 2  m 4  f  4  p 3  r 3

(7)

All coefficients in Eq.(7) are material constants, which have physical meanings and are related to certain material properties, but there are not enough experiments to fit or otherwise estimate all these material constants.. Hence, simplifications are necessary. Based on the above information, we assume that the free energy W1 is a mathematical transformation from W in the Eq. (8) (Khachaturyan, 2013). First, we consider the non-convex free energy

W  ,   , where  denotes the twist (per unit length) and  denotes the curvature (per unit length); then W can be expressed as

W  a 2  b 4  d 2

(8)

 and  serve to define the minimum values of the energy function in the normal and

semi-coiled

forms respectively. Accordingly, the energy function has local minima at points

(  1 = – 2.0rad/μm, 1 = 1.2rad/μm) and (  2 = 2.0rad/μm  2 = 2.4rad/μm), as noted previously.

dW   ,   dW   ,    1, 1  0,  1, 1  0 d d d 2W   ,   d 2W   ,   d 2W   ,   A  1, 1 , B   1, 1 , C   1, 1 d 2 d  d 2

(9)

  AC  B 2  0, A  0. From the equation (9) we obtain a  8.2; b; and d  0 ; experimental data are required to 7

determine the exact values of a, b and d. Secondly, according to Flynn (Flynn and Ma, 2004) and Hoshikawa (Hoshikawa and Kamiya, 1985), in the normal form, the dimensionless twist/bend ratio (EI/GJ) is in the range of 0.005 to 5.0, and so, our result W /  W / 

 0.4 is within that range. Finally, from the analysis of the experimental work (  N , N )

of Hotani (Hotani, 1982a), we see that the torque and bending for the viscous loading is on the order of 10−19 Nm, from these values we estimate that b is about 10−25 Nm2, and d is equal to a. So, we set a  d  8.2 1025 , b  11025 . The constants are comparable with the experimental results of Fujime (Fujime et al., 1972): in that study it was determined through quasi-elastic scattering of laser light that the flexural rigidity of the straight bacterial filament is about 10 15dyncm2 1024 Nm2,

. In the

model used in the current study the flexural rigidity is about 2 10-25 Nm2, the torsional rigidity is about 3.3 10-25 Nm2. Also based on Hotani’s experimental results, the flexural rigidity along the e2 direction is estimated to be order of 1023 Nm2. We assume that the mathematical transformation of the geometry includes translation and rotation of the original energy function, as shown in Figure 5(a). We assume that the transformation of the geometry involves translating along the x direction by about Δx, translating along the y direction by about Δy and rotating around z by about θ, and obtain the following relationship between the two coordinates, as shown in Figure 5(b).

   x'   x    '   y   y        '    T ; xy     R;    T ; x, y        z    z       1    1        

(10)

By going through the above mathematical transformation from Equation (10), we finally 8

get the elastic energy density of the filament as equation (11). The energy map is plotted in Fig. 6, showing two energy wells corresponding to the two stable phases of filament, namely the left handed normal and the right handed semi-coiled. For mathematical convenience, we shift the twist-rates of the two phases from (2.0 rad/μm, 2.0 rad/μm) to (0 rad/μm, 4 rad/μm and the curvatures from (1.2 rad/μm, 2.4 rad/μm) to (0 rad/μm, 1.2 rad/μm). W  0.3a  0.1b  2.9d  a  0.6b  d   0.9a 2  1.6b 2  0.1d  2  1.9b 3 0.8b 4  0.3a  0.2b  3.3d  0.6a  b  0.6d   1.8b 2  b 3 a  0.2b  d  0.6b  0.5b   0.1b  0.1b 2

2

2

2

2

2

3

(11)

3

An important simplification to the energy expression given above results in the case where the rod is approximately straight: Once a choice is made for a given configuration, such as the stress-free state, that choice also applies to all other configurations. In the classical rod theory, it is natural to align e1and e2 with the principal axes of the cross-section. For example, if the rod aligns along the z-axis when it is stress free, then we use the coordinate system

 x, y,z

e1 ,e2   x, y

is natural. Here

where the un-deformed filament is aligned with the

z-axis, as put forward by Landau and Lifshitz (Nowacki, 1986). The position of the filament is now written as the triad

 X  z  , Y  z  , z We also let θ(z) be the angle through which the

material frame has been twisted as a function of z; then the filament configuration is completely determined by

 X  z  , Y  z  ,   z 

Therefore, to the same order in Xz, Yz we

simply insert linearized approximations for  , and Ω2 into the energy. The material frame now rotates with θ(z) and so the linear strain approximations are rotated appropriately per equation (12)

9

2 X  2Y  sin  z 2 z 2 2 X  2Y  2   sin  2  cos  2 z z    3  z

  cos 

(12)

2.3 2D non-local energy

A general method to analyze the interfacial properties was proposed by Cahn and Hilliard (Cahn and Hilliard, 1958), who suggested that the free energy of each point of a non-uniform system is not only related to the component of this point, but also depends on the distribution of the component in the small area surrounding the point. Using Taylor series, the free energy of each point can be expressed in two parts, one part being the local energy which depends on the local value of the component at the point, and the other part being the gradient energy which represents the change in the value of the component near the point (Mindlin, 1965; Mindlin and Eshel, 1968; Triantafyllidis and Aifantis, 1986; Triantafyllidis and Bardenhagen, 1993). Here, in order to describe the interface property, we introduce the strain gradient energy as an additional term in the total system free energy.

At the interface, the corresponding transformation strains are interpolated by a continuously varying elastic strain field that depends only on the position variable ―z in the direction of the filament, by analogy to figure7 for the one-dimensional ―φ-4 model (Barsch and Krumhansl, 1992; Olson and Owen, 1992). The variation of the two-component order parameter

 , 

as a function of position s may then be visualized qualitatively in terms of

a curve     s  ,    s  that moves in a three-dimensional  , ,s  coordinate system (s moves along the filament) from tetragonal energy minimum I to the minimum II. The 10

non-local strain gradient energy density can be expressed as

G  g1 (

 2    )  g 2 ( ) 2  g3 s s s s

(13)

with the derivatives expanded in terms of the material coordinates as   cos  X zzz  sin  z X zz  sin  Yzzz  cos  zYzz z    zz z

(14)

where g1 , g2 , g3 contain a material intrinsic length characterizing the interface thickness. 3. FEM methodology The procedures for implementing the non-convex and non-local model in FEM simulation are as follows, with the assumption that the filament is inextensible. In order to integrate the free energy along the filament (x coordinate) in Eq. (1), the filament is discretized into elements. The displacement field (distribution of the rotation angle and bending deflection) in each element is represented by the nodal displacements and proper shape functions. In the present FEM model, the nodal displacements are the rotation angle and the bending deflection; the derivatives of the displacement are the twist and the curvature. As shown in Eqs. (1) and (14), the energy of the system includes twist, twist gradient, curvature and curvature gradient. So elements with C-2 continuity are used in the current model. In the FEM simulation shown in Fig. 8, the filament (of length 20 μm) is divided into n-1 elements with n nodes. Each element has two nodes and each node has five nodal degrees of freedom: rotation angle  , twist  ,x , bending deflection u, curvature u,x, curvature gradient u,xx. Then the displacement field  in an element can be calculated as 11

i     i,x  ui    ui , x  u   i , xx     h1 h 2 h3 h4 h5 h6 h7 h8 h9 h10     i 1  i 1, x    ui 1  u   i 1, x  ui 1, xx   

(15)

where the shape functions are

h1 = h2 = h3 = h4 = h5 =

1 (1-m)3 (1+m)(5+3m), 32 h2 h7 = (1-m)3 (1+m) 2 , 64 1 h8 = (1+m)3 (8-9m+3m 2 ), 16 h h9 = (1+m)3 (-1+m)(5-3m), 32 h2 h10 = (1+m)3 (1-m) 2 . 64

1 (1-m) 2 (2+m), 4 h (1-m) 2 (1-m), 8 1 (1+m) 2 (2-m), 4 h (-1+m 2 )(1+m), 8 1 (1-m)3 (8+9m+3m 2 ), 16

h6 =

Here, m is the local coordinate of an element; and the local coordinates of the two nodes in one element are -1 and 1, respectively (Fig. 6). With the displacement field (rotation angle field and bending deflection field) in Eq. (15), the twist  , twist gradient  ' , curvature

 and curvature gradient  ' can be expressed as  h1 h2 x

   x

  2 h1  2 h2 2 x 2  x

 '

h6 x

i    h7  i , x    x  i 1  i 1, x    i     2 h7  i , x    x 2  i 1  i 1, x   

 2 h6 x 2

12

(16)

(17)

 h3

   x

  2 h3  x

'

h4 x

 2 h4 x

h5 x

h8 x

 2 h5 x

 2 h8 x

h9 x

 2 h9 x

ui  u   i,x  h10  ui , xx    x  ui 1  ui 1, x    ui 1, xx 

(18)

ui  u   i,x   2 h10  ui , xx    x  ui 1  ui 1, x    ui 1, xx 

(19)

Substituting Eq. (16), (17), (18), (19), (11), (12) and (13) and (11)-(14) into the energy variation Eq. (3) and assembling all the elements, we can obtain (Sun and He, 2008; He and Sun, 2010) the equilibrium equations as ~

~

K q  0

(20)

~

~

where q is a vector containing all the nodal freedoms of the system; K is the stiffness. Since the material constitutive energy is nonlinear (i.e. non-convex elastic energy in Eq. (15), ~

K is not a constant matrix, but depends on the state of the material (i.e., depends on the ~

nodal freedoms q ). Therefore, the equilibrium equations (Eq. (20)) are nonlinear.

The

Newton-Raphson iteration method (He and Sun, 2010) is used to solve such nonlinear equations of the form       0,

(21)

where  is a nonlinear function of the unknown q while  is a constant. If an approximation (initial) value of the unknown is given, q(i), an improved solution, q(i+1), can be obtained through a Taylor expansion of Eq. (22) for the first order approximation: 13

 (q(i 1) )   (q(i ) ) 

 (q(i ) ) q

.q(i )  0

 1 )  (q(i ) ) q  (   ) 1  ( )  (q(i ) ) q   ( ) 1 (q(i ) ) q  q(i 1)  q(i )  q(i ) .  (

(22)

This means that the value in the (i+1) th iteration step, q(i+1), can be obtained by that in the ith iteration step, q(i), through Eq. (23). The iteration process is repeated until the following convergence criterion is met

{[q(i 1) ]}2  a small value.

(23)

The displacements δ (rotation angle and bending deflection) at the two ends of the 2D system (n nodes) are specified as 1  0 , and  n   ), where  is the specified proportional loading between rotation angle and bending deflection with 0.4. The small circle on the energy curve in Fig. 3 is the spinodal point of the non-convex energy. When a perfect filament is loaded up to this point, the FEM simulation encounters an ill-conditioned stiffness and the iterations fail to converge. An artificial defect placement is a general way to avoid this bifurcation point. In our current simulations, several selected elements are hence treated as defects by reducing their cross-section area by a certain percentage so that new phase can be nucleated there. We set eight out of all eighty elements as defects at the fixed end, by reducing their cross-section area to 10% of that of the other elements, thus ensuring that the new phase will be nucleated at these defects. 4. Results and Discussion

14

We simulate the filament phase transition from the (initial) N phase to the SC phase, with the filament fixed at one end, and the proportional loading between rotation angle and bending deflection set at 0.57 the other end. The number of loading steps is 500 for total torsion angle 25 rad, and the total bending deflection is 16.2 μm. The torque-twist curve and bending-curvature curve of a filament (length 20 μm) at the loading end are shown in Fig. 9(a) and 9(b). After the first elastic deformation of the N phase, both the torque and bending drop to static values, which correspond to the SC phase nucleation and the new-phase propagation. During the growth of the SC phase, both torque and bending remain static, which means that the two phases coexist, as shown in Fig. 9 (e) and (f). When the SC phase grows to occupy the whole filament, the torque and bending experience a sudden change due to interface annihilation. Further loading leads to elastic deformation of the SC phase. Corresponding to the response of torque Fig. 9(a) and bending 9(b), the twist and curvature profiles along the filament of the phase nucleation and growth are shown in Fig. 9(e) and 9(f). After the initial elastic loading (loadings 1-3), an SC phase nucleates from the N phase (loading 4). Then the SC phase grows and coexists with the N phase. The equilibrium twists of these two phases are near 0 rad/μm and 4 rad/μm, while their equilibrium curvatures are near 0 rad/μm and 1.2 rad/μm respectively. There is a transition zone or interface between the two phases, of thickness around 1 μm. During the coexistence of the two phases, the SC phase grows via the interface propagation when the proportional loading of rotation angle and bending deflection increases (loadings 4-6 in Fig. 9(e) and 9(f)). Loading-unloading cycles of rotation angle and bending deflection are shown in Fig. 9(c) and 9(d). Both torque and bending drop down in the loading curve (signifying phase

15

nucleation, with the arrow down) and jump up in the unloading curve (signifying phase annihilation, with the arrow up). It is noticeable that a loop is formed by the loading and unloading curves, which means that some energy is dissipated as elastic vibration and heat during the unstable processes (phase nucleation and annihilation). The twist and curvature profiles corresponding to the response of the proportional unloading of rotation angle and bending deflection are shown in Fig. 9(g) and 9(h). With the proportional unloading, the SC phase shrinks via interface propagation (unloading steps 1-3). When the phase shrinks to a critical small size (about twice the interface thickness), the phase is annihilated suddenly (unloading steps 4-6). Then the filament returns to the initial N phase. If we decrease the gradient energy by reducing g from 1 μm to 0.08 μm, the calculated interface thickness is comparable with what is experimentally observed (0.2-0.08 μm (Coombs, 2001)); we obtain the corresponding twist profiles as Fig. 9(e) (g = 1 μm) and Fig. 6(a) (g = 0.08 μm), and the curvature profiles as Fig. 9(f) (g = 1 μm) and Fig. 10(b) (g = 0.08 μm), respectively. We can see that the interface thickness decreases with decreasing g. For simplification, we set g=g1=g2=g3. The drop in the interface thickness means that the interface occupies less volume and has less interfacial energy. This must affect the processes of phase nucleation and phase annihilation. In other words, the changes in the interface properties (governed by g in the current model) affect the phase transition. The phase nucleation location is not merely the initial point. In order to further study the multiple phases nucleate in the filament, we use the 1D model schematically to simulate the mechanical behavior of the filament. In our current simulations, with the total number of elements set at 80, we set two locations of defects from elements 28 to 35 and elements 48

16

to55, by reducing the cross-section area to a percentage of 10%, so that a new phase can be nucleated at the defects. The filament is fixed at node 1, while the rotation angle loading is added at node 81. The torque-twist curve of a filament (of length 20 μm) is shown in Fig. 11(a). After the first elastic deformation of the N phase, there is a torque drop in the curve which corresponds to the nucleation of SC phases at the defect locations. During the growth of the SC phases, the torque is zero, which means that the two phases coexist on the Maxwell line (Wang and Sun, 2010a). When neighboring interfaces of two SC phases merge, the torque experiences a sudden drop due to interface annihilation. Corresponding to the response of torque Fig. 11(a), the twist profile of the phase nucleation and growth is shown in Fig. 11(b). After the initial elastic loading (loadings 1-2), the SC phase nucleates simultaneously from the N phase in two locations (loading 3). Then it grows and coexists with the N phase. The equilibrium twists of these two phases are near 0 rad/μm and 4 rad/μm, respectively. There is a transition zone (interface) between the two phases, of thickness around 1 μm. During the coexistence of the two phases, two SC phases propagate to the two ends separately (loadings 3-4), then one end of one SC phase merges with the neighboring end of another SC phase (loading 4). During the growth of the SC phase, both torque and bending remain static, which means that the two phases coexist. Thus the number of interfaces decreases from 4 to 2 due to the interface merger and annihilation. Finally the merged SC phase propagates along the filament when the rotation angle loading increases (loading 5-6). A loading-unloading cycle is shown in Fig. 11(a). The torque drops down (phase

17

nucleation, with the arrow down) in the loading curve while the torque jumps up (phase annihilation, with the arrow up) in the unloading curve at the smaller twist; a loop is formed by the loading and unloading curves, which again means that some energy is dissipated as elastic vibrations and heat during the unstable processes (phase nucleation and annihilation). At the lager twist, between the torque drop in the loading curve and torque jump in the unloading curve, there is a smaller loop, which is shown more clearly in the inset of Fig. 11a. This loop signifies a reduction in the system energy that is mainly caused by the decrease in the interfacial energy of the domain during the merger (annihilation) of neighboring interfaces when the two SC phases grow into one; the energy is lost to elastic vibrations and heat. A similar phenomenon of energy dissipation caused by instability due to phase transition and phase annihilation was found in shape memory alloys (He and Sun, 2009a).

Conclusion In summary, we introduce a 2D non-convex non-local viscoelasticity model and implement it in FEM simulations to describe the flagellar phase transition between the N phase and the SC phase. The total free energy in our model comprises three parts: The non-convex elastic energy describes the material instability involved in the phase transition; the non-local gradient energy governs the thickness and the energy of the interface (transition zone) between the two coexisting phases; further, viscous dissipation is included to account for the relaxation process associated with non-equilibrium conditions. In all, this model and the FEM code provide a basic tool to investigate the behaviors of flagella under various loading conditions. In actual (as opposed to simulated) experiments, the flagella phase transition can be very complicated, e.g. cyclic in a flowing fluid, with a deformation that

18

includes both twist and curvature. So, in order to achieve a better understanding of the flagellar behavior, we extend our previous work to 2D by considering two order parameters in free energy.

The main conclusions of the paper are as follows: (1) A non-convex nonlocal elasticity model and the associated FEM code in two dimensions are developed to simulate the transformation process of filament phase nucleation, propagation and annihilation. The results of the 2D model are more consistent with the experimentally observed phenomena (torque and bending accumulation which correspond to phase nucleation, torque and bending drop which correspond to phase growth, and forward and backward transitions) (Wang and Sun, 2010a). (2) Transition instability, i.e. the switching between different equilibrium configurations, causes energy dissipation and a change in the number of interfaces. The loop in the loading and unloading curve indicates energy dissipation. Moreover, there is load drop when the two interfaces merge, which means that a reduction in the number of interfaces reduces the interface energy. (3) The nonlocal interface energy can describe the properties of interface between the lefthanded normal form and the right-hand semi-coiled form. The interfacial properties affect both the twist and curvature profiles as well as the torque and bending responses during phase nucleation and phase annihilation. The interface between the two phases is actually a transition zone, and interface thickness is governed by the coefficients of the gradient of order parameters.

19

Acknowledgements The authors would like to thank Harvard University for their support. The National Natural Science Foundation of China (11772047), and Key international collaborating Project from National Natural Science Foundation of China (11620101001).

Conflict of interest The authors declare that they have no conflict of interest.

References Asakura, S. 1969. Polymerization of flagellin and polymorphism of flagella. Advances in biophysics 1, 99-155. Asakura, S., and Iino, T. 1972. Polymorphism of Salmonella flagella as investigated by means of in vitro copolymerization of flagellins derived from various strains. Journal of molecular biology 64, 251IN17257-17256IN25268. Barsch, G., and Krumhansl, J. 1984. Twin boundaries in ferroelastic media without interface dislocations. Physical Review Letters 53, 1069. Barsch, G., and Krumhansl, J. (1992) Topological soliton models for interfaces in proper ferroelastic martensites. Proceedings of the international Conference on martensitic transformations. Cahn, J.W., and Hilliard, J.E. 1958. Free energy of a nonuniform system. I. Interfacial free energy. The Journal of chemical physics 28, 258-267. Calladine, C. 1976. Design requirements for the construction of bacterial flagella. Journal of theoretical biology 57, 469-489. Calldine, C. 1978. Change of waveform in bacterial flagella: the role of mechanics at the molecular level. Journal of Molecular Biology 118, 457-479.

20

Coombs, D. 2001. Dynamics of travelling helicity fronts in bacterial flagella. Coombs, D., Huber, G., Kessler, J.O., and Goldstein, R.E. 2002. Periodic chirality transformations propagating on bacterial flagella. Physical review letters 89, 118102. Curnoe, S., and Jacobs, A. 2001. Time evolution of tetragonal-orthorhombic ferroelastics. Physical Review B 64, 064101. Darnton, N.C., and Berg, H.C. 2007. Force-extension measurements on bacterial flagella: triggering polymorphic transformations. Biophysical journal 92, 2230-2236. Darnton, N.C., Turner, L., Rojevsky, S., and Berg, H.C. 2007. On torque and tumbling in swimming Escherichia coli. Journal of bacteriology 189, 1756-1764. Falk, F. 1983. Ginzburg-Landau theory of static domain walls in shape-memory alloys. Zeitschrift für Physik B Condensed Matter 51, 177-185. Flynn, T.C., and Ma, J. 2004. Theoretical analysis of twist/bend ratio and mechanical moduli of bacterial flagellar hook and filament. Biophysical journal 86, 3204-3210. Fujime, S., Maruyama, M., and Asakura, S. 1972. Flexural rigidity of bacterial flagella studied by quasielastic scattering of laser light. Journal of molecular biology 68, 347IN3355-3354IN4359. Goldstein, R.E., Goriely, A., Huber, G., and Wolgemuth, C.W. 2000. Bistable helices. Physical Review Letters 84, 1631. He, Y., and Sun, Q. 2009a. Effects of structural and material length scales on stress-induced martensite macro-domain patterns in tube configurations. International Journal of Solids and Structures 46, 3045-3060. He, Y., and Sun, Q. 2009b. Scaling relationship on macroscopic helical domains in NiTi tubes. International Journal of Solids and Structures 46, 4242-4251.

21

He, Y., and Sun, Q. 2010. Macroscopic equilibrium domain structure and geometric compatibility in elastic phase transition of thin plates. International Journal of Mechanical Sciences 52, 198-211. Hoshikawa, H., and Kamiya, R. 1985. Elastic properties of bacterial flagellar filaments: II. Determination of the modulus of rigidity. Biophysical chemistry 22, 159-166. Hotani, H. 1976. Light microscope study of mixed helices in reconstituted Salmonella flagella. Journal of molecular biology 106, 151-166. Hotani, H. 1982a. Micro-video study of moving bacterial flagellar filaments. J. Mol. Biol 156, 791-806. Hotani, H. 1982b. Micro-video study of moving bacterial flagellar filaments: III. Cyclic transformation induced by mechanical force. Journal of molecular biology 156, 791-806. Kamiya, R., and Asakura, S. 1976a. Flagellar transformations at alkaline pH. Journal of molecular biology 108, 513-518. Kamiya, R., and Asakura, S. 1976b. Helical transformations of Salmonella flagella in vitro. Journal of molecular biology 106, 167-186. Khachaturyan, A.G. 2013 Theory of structural transformations in solids. Courier Corporation. Lazarus, D., and Jahn, R. 1998. Using the Ehrenberg collection. Diatom Research 13, 273-291. Mindlin, R., and Eshel, N. 1968. On first strain-gradient theories in linear elasticity. International Journal of Solids and Structures 4, 109-124. Mindlin, R.D. 1965. Second gradient of strain and surface-tension in linear elasticity. International Journal of Solids and Structures 1, 417-438. Nowacki, W. 1986. Theory of asymmetric elasticity. Pergamon Press, Headington Hill Hall, Oxford OX 3 0 BW, UK, 1986. Olson, G.B., and Owen, W.S. 1992 Martensite: a tribute to Morris Cohen. Asm Intl.

22

Reid, A., and Gooding, R. 1994. Hydrodynamic description of elastic solids with open boundary conditions undergoing a phase transition. Physical Review B 50, 3588. Srigiriraju, S.V., and Powers, T.R. 2005. Continuum model for polymorphism of bacterial flagella. Physical review letters 94, 248101. Srigiriraju, S.V., and Powers, T.R. 2006. Model for polymorphic transitions in bacterial flagella. Physical Review E 73, 011902. Sun, Q., and He, Y. 2008. A multiscale continuum model of the grain-size dependence of the stress hysteresis in shape memory alloy polycrystals. International Journal of Solids and Structures 45, 3868-3896. Triantafyllidis, N., and Aifantis, E.C. 1986. A gradient approach to localization of deformation. I. Hyperelastic materials. Journal of Elasticity 16, 225-237. Triantafyllidis, N., and Bardenhagen, S. 1993. On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models. Journal of Elasticity 33, 259-293. Turner, L., Ryu, W.S., and Berg, H.C. 2000. Real-time imaging of fluorescent flagellar filaments. Journal of bacteriology 182, 2793-2801. Wada, H., and Netz, R.R. 2008. Discrete elastic model for stretching-induced flagellar polymorphs. EPL (Europhysics Letters) 82, 28001. Wang, X.-L., and Sun, Q.-P. 2010a. Mechanical analysis of phase transition experiments of the bacterial flagellar filament. Acta Mechanica Sinica 26, 777-785. Wang, X. 2006. Mechanical analysis and free energy construction of phase transition in bacterial flagellar filaments.

23

Wang, X., He, Y., and Sun, Q. 2011. Simulation of bacterial flagellar phase transition by non-convex and non-local continuum modeling. Theoretical and Applied Mechanics Letters 1, 044001. Wang, X., and Sun, Q. 2010b. Phase transition in biological systems. Adv. Mech 40, 41-56.

24

FIGURE LEGENDS Fig. 1 (a) Illustration of movements of Escherichia coli, when its motor rotates CCW, flagella are in a bundle and the cell swims; when the motor reverses direction, flagella fall apart and the cell tumbles. When the motor reverses again, flagella re-bundle and the cell swims in another direction. (b) The experiment of Hotani (Hotani, 1982) captures external viscous fluid force induced filament phase transition; namely, the cyclic transformation between normal and semi-coiled forms. A flagellar filament is attached to a glass surface at its top portion and is subjected to a flow of viscous fluid from top to bottom. Successive frames (from left to right) of the cinemicrograph are taken at intervals of 5/9s. The bar represents are 10 μm. (c) A sketch of figure 1(b) illustrating the cyclic transformation between normal and semi-coiled form under external viscous fluid force. Figure 1 is reprinted with permission from Ref.(Wang and Sun, 2010a)

Fig. 2 (a) All component forces generated on the helical filament. Bending moment M1, M2; torque T; shear force, Q1, Q2; tension, F. (b), (c), (d) or (e) Torque T (M1, M2 or F) (Nm or N, respectively) at the filament cross section, The black lines with black dots are T (M1, M2 or F) at stuck end respectively; the color lines with black dots are T (M1, M2 or F) at the forward transition (normal to semi-coiled) points; the black lines with color dots are T at the reverse transition (semi-coiled to normal) points. (Wang and Sun, 2010a)

Fig. 3. a). shows a filament with a stripe in a straight, untwisted state. b). shows the curvature perpendicular to the stripe, that is around the local e2 axis. c). shows the curvature around the local e1 axis; d). shows twist around the tangent e3 axis.

Fig. 4 Space natural frame: tan is along the tangent direction, n is along the curvature direction called normal direction and b is along the bi-normal direction. t indicates time.

Fig. 5 (a)The different order of rotation and translation results in different transforming picture; (b) The mathematical expression for different kinds of graph transformation. The above one is for the order of rotation and translation; the lower one is for the order of translation and rotation.

Fig. 6 Double-well free energy as a function of twist and curvature. The two wells represent two stable phases: normal and semi-coiled. The black circle represents the spinodal point

Fig.7 Schematic diagram showing (dashed line) the variation of an order parameter φ as a function of position s for the “φ-4” model. This topological, localized kink connects phase I and II. V(φ) is the 2 4 non-convex potential    a / 2     b / 4    .

Fig. 8 Nodal freedoms and elements used in the FEM simulation. The filament is divided into n-1 elements; each element has five nodal freedoms:

i ,i , x , ui , ui , x , ui , xx ; h is the element length.

Fig. 9 (a) Torque-twist curve of a filament under proportional loading of rotation angle and bending 25

deflection; (b) Bending-curvature curve of a filament under proportional loading; (c) Torque-twist curve of proportional loading and unloading of rotation angle and bending deflection; (d) Bending-curvature curve of proportional loading and unloading; (e) Twist profile (g = 1 μm) of a filament under proportional loading; (f) Curvature profile (g = 1 μm) of a filament under proportional loading; (g) Twist profile of a filament under proportional unloading; (h) Curvature profile of a filament under proportional unloading.

Fig. 10 (a) Twist profile (g = 0.08 μm) under proportional loading; (b) Curvature profile (g = 0.08 μm) under proportional loading.

Fig. 11 (a) Simulated torque-twist curve of a filament (two defects) under loading and unloading of rotation angle; (b) Simulated twist profile of a filament (two defects) under loading of rotation angle.

26

a

b

Attached end

Normal

Semi-coiled

Fluid direction

10 μm c

Fluid direction

Semi-coiled

Time interval 5/9s

Rotation direction

Fig.1

27

Normal form

Stuck end

a

b

c

d

e

f

Fig.2

28

a b

c

d

Fig.3

Fig.4 a

b

Fig.5 29

Fig.6

Fig.7

30

Fig.8

31

Fig.9 32

Fig.10

Fig.11

33