Shear-induced microstructures and dynamics processes of phospholipid cylinders in solutions

Shear-induced microstructures and their corresponding dynamic processes are investigated for phospholipid cylinders in aqueous solution by dissipative particle dynamic simulation. Various phospholipid cylinders with cross-sections, which are formed under shear-free flow, are selected to examine the effects of shear flow on their structures and dynamic processes. Shear flow induces the transition from cylinders into vesicles at weak rate and the transition into vesicle–lamella mixtures with increased shear rate and lamella structures at the strong shear rate. Then, the average radius of gyration and shape factors of the polymer chains in the dynamic processes are discussed in detail. Results show that shear flow causes the structure of the polymer chains to be elongated along the shear direction, and the configuration of the polymer chain can be rapidly transformed into an ellipsoid structure under strong shear.

the cylindrical network structure can be observed (Fig. S1 in the Supplementary Materials), and this condition is independent of chain lengths in the concentrated solutions 14 .
Applying shear flow provides an external force that can induce the phase transitions of phospholipid polymers in solutions instead of varying the inner polymer parameters, such as chain lengths; as a result, corresponding dynamic processes are introduced into the microstructures 15,16 . Phospholipid membranes, a basic bio-structure in cells exposed to hydrodynamic fluid flows, under shear flows have been increasingly investigated, such as hemoglobin flowing in the blood. Shear flows not only redistribute the molecular locations in membranes but also induce the structural changes of membranes and the corresponding dynamic processes caused by hydrodynamic forces [17][18][19][20][21][22] . For example, external liquid flows result in hydrodynamic forces on molecules on the phospholipid bilayer that causes them to move in the direction of the shear flow and leads to the local concentration of molecules change. Quantitative studies on the magnitude of shear forces have been carried out 21 . MD simulation predicted that shear force leads to the rotation and re-alignment of molecules in phospholipid membranes; these phenomena result in structural undulations that propagate in the perpendicular direction when the shear rate exceeds a critical value 19 . However, layer sliding can occur instead of undulation in phospholipid bilayers under shear flow, and the intermonolayer slip depends on the velocity profile 20 . Dynamic processes, such as changes in system energy and changes in the radius gyration of structures, have also been investigated. For example, the system energy and radius of gyration have been calculated by MD simulation for phospholipid membranes; the results showed that system energy tends to be stable, and gyration radius increases when the shear increases 14 . Shear flows also induce various microstructures, not just lamellar structures, for block copolymers in solutions 23 . In particular, shear flows can induce structural changes in phospholipid vesicles whose shapes are similar to red blood cells, which become exposed to hydrodynamic fluids and cause structural changes [24][25][26][27][28][29][30][31][32][33] . Additionally, shear flows can induce shape changes in entire phospholipid vesicles along with domain diffusions and round domains 24 ; the deformation of phospholipid vesicles into steady ellipsoidal shapes in constant orientations under simple linear shear flow has been observed by phase contrast microscopy 28 .
Many experiments and simulations have been performed to determine the effects of shears on phospholipid membranes and vesicles. However, the shear effects on phospholipid cylinders still need to be explored. In the current study, we investigate the shear-induced microstructures and the dynamic processes of phospholipid cylinders in solutions by using the DPD and CG model, which is useful in simulating complex fluids and soft matters [34][35][36][37] . We focus on the influence of shear flow on phospholipid cylinders with various cross-sections in different concentrations. By using a variety of initial state inputs, the energy of the formed equilibrium structure can be compared, and the simulation method of the stable structure can be selected. Different shear rates are considered to examine the effects on the microstructures and their corresponding dynamic processes. The rest of the paper is structured as follows. Section II describes the model and the method. Section III presents the results and discussion. Section IV presents the summary.
Model and Method phospholipid model. We coarse-grain a small group of atoms into a single bead in the current phospholipid model. The phospholipid molecule is modelled with two-tail linear chains and one head linear chain, which has been extensively used in the previous studies, as shown in Fig. 1a 38 . In the phospholipid model, the hydrophilic heads (H) and hydrophobic tails (T) are represented by red and yellow beads, respectively. In order to connect two consecutive beads, we used an additional elastic harmonic force in the model. And this elastic harmonic force can be expressed as follows: www.nature.com/scientificreports www.nature.com/scientificreports/k where the ith and jth beads are two consecutive beads, and k s represents the spring constant and r s is the equilibrium bond length between the two beads, whereas r ij denotes the distance between ith and jth beads, with the unit vector r ij . We take k s = 100.0 and r s = 0.7 r c , which are similar to those in the previous studies 39,40 . Here, we set r c is the cutoff radius in the DPD simulations. An additional force caused by the harmonic constraint is applied onto two the consecutive bonds to achieve the bending resistance of the chains and is given as follows: where k θ , θ and θ 0 are the bending constant, inclination angle and equilibrium angle, respectively. We take θ 0 = π for three consecutive T beads and θ 0 = 2/3π for the three consecutive beads at the connective point between the T and H, as shown in Fig. 1a, which are similar to those of the model used in our previous works 13,41 . In this study, we fix N HB = 3 and N TB = 4 for H and T, respectively. The model phospholipids are integrated into the aqueous solution, in which the water molecules are also coarse-grained. Then, we define the phospholipid concentration as follows: Shear flow. In linear response theory, a flux relates to a thermodynamic force or field. In the isotropic fluid, the flux has a linear relationship with the corresponding field with a transport coefficient. In non-equilibrium MD, the field is applied, whilst the flux is measured 42,43 . By contrast, in reverse non-equilibrium MD, the momentum flux is introduced, whilst the force or field is measured 18,44,45 . We have successfully applied this reverse non-equilibrium approach in the DPD simulations 14 . Here, we only briefly discuss this approach as follows. The momentum flux is introduced into the system in an unphysical way. Specifically, the periodic simulation box with box length L z is divided into several slabs in the z-direction. The beads inside the slabs at z = 0 and z = L z are propelled in two opposite x-directions, respectively, as shown in Fig. 1b. The particles can then obtain the largest momentum in the opposite directions in the two slabs. In proving that such shearing is reasonable, the velocity <v x > distribution in this method under the three shear rates is used to roughly satisfy the linear relationship, as shown in Fig. 2 44 . Then, the momentum component p x is exchanged by executing the flip algorithm at regular time intervals. Thus, the momentum flux j z (p x ) can be calculated by where t is the time length of the two swaps, and factor 2 arises because of the periodicity of the simulation box with L x and L y in the x-and y-directions, respectively. During simulation, we perform the velocity swap every W time steps to satisfy t = WΔt. We can adjust the time lapse WΔt between the two velocity swaps to control the momentum flux j z (p x ). In this study, we apply the steady shear flows, in which the velocity distributes linearly in the z-direction. In the steady state, the rate of momentum transferred by the momentum swaps is equal to the momentum flowing back through the fluid by friction. In this manner, the unphysical momentum swap conserves the linear momentum and the kinetic energy of the system as a whole because bead positions are unchanged. Moreover, shear rate can be calculated from with the collinear momentum flux in this method. On this basis, we can obtain the shear rates / 1 γ τ − = 0.02, 0.04, 0.06 by setting W at 5, 2 and 1, respectively.
DpD method. In order to simulate the hydrodynamic behaviour of complex fluids and soft matters, we often use the DPD method. In DPD simulations, several pairs of forces are introduced to utilise beads that are coarse-grained from a small group of atoms. The main formulations about the DPD simulation are similar to those used in complex fluids and soft matters 35,39,[46][47][48] . During DPD simulation there are three forces acting on the ith bead: a conservation force (F ij C ) from a potential, a dissipative force (F ij D ) that attempts to reduce radial velocity differences between the particles and a random force (F ij R ) that represents stochastic impulse. We can expressed the total forces on the ith bead as follow: where a ij is the maximum repulsive force between the ith and jth beads, and r ij is the distance between the ith and jth beads and v ij denote the relative velocity between the ith and jth beads, with the unit vector (r iĵ ). We take a ij = 200 for the hydrophilic type of beads and a ij = 50 for the hydrophobic types. The validity of fluctuation-dissipation theorem requires two parameters, namely,  γ and σ, to be linked by where k B is the Boltzmann constant and T is the system temperature 49 . We used σ = 3.0 in simulation which is usually the standard valued before 39 . ζ ij represents a Gaussian distribution of zero mean and unit variance, and the weight function w (r ij ) can be expressed as where r c is the cutoff radius.
Simulation parameters. The length, time and energy in DPD simulation can be normalised by several parameters. In particular, energy can be normalised by k B T, whilst length is normalised by the cutoff radius r c in accordance with r c = (ρV b ) 1/3 , where V b is the volume of one DPD bead, and ρ is the bead density. Here, we set ρ = 3 in the simulation, and the DPD bead usually has an assumed volume of Vb = 0.03 nm 3 , which leads to a value of r c = 1.0 nm 50,51 . Time unit τ is defined as τ = r m k T / c B , where m is the bead mass. The time unit can be deduced by considering the in-plane diffusion constant of the lipid in the experiment 52,53 . The movement of the DPD beads is driven by DPD forces based on Newton equations, which are integrated by a modified version of the velocity Verlet algorithm with the time step of Δt = 0.01τ 47,54 . All simulations are performed on an NVT ensemble by using LAMMPS. The simulation is performed in a box set to V = L × L × L, we set all three directions with periodic boundary conditions 34 . Additional simulations are carried out by varying the box sizes to avoid the finite size effect, and the box sizes are optimised to L = 30r c 55 . After inputting different initial states, the self-assembly is performed, the energy value in the equilibrium state is compared and the structure with the smaller energy value is selected as the final stable structure. According to the comparison of simulated energy, the final stable state can be achieved after 200000 DPD time steps, which is similar to those in the previous works 56,57 .

Results and Discussion
We mainly focused on the microstructures of phospholipid porous cylinders with chain lengths of N HB = 3 and We observe porous cylinder structures with three small cylinders in a triangular arrangement at φ P = 0.35, with four small cylinders in a quadrangular arrangement at φ P = 0.40, with five small cylinders in pentagonal arrangement at φ P = 0.45 and with six small cylinders in hexagonal arrangement at φ P = 0.50, as shown in Fig. 3. In other words, the phospholipid concentration causes the small cylinders to be packed inside the large cylinder. This porous cylinder differs from the cylinder observed in the previous experiment and simulation, in which the cylinders are single or hexagonally arranged without porous cross-sections in the solutions 3,8 . Here, we observe phase transitions between the porous cylinders caused by the phospholipid concentrations.
The microstructures have weak shear rates of 0.02, as shown in Fig. 4. Considering that shear flow is applied perpendicularly onto the axes of the porous cylinder, structural transitions occur even under weak shear flows. In (2019) 9:15393 | https://doi.org/10.1038/s41598-019-51933-z www.nature.com/scientificreports www.nature.com/scientificreports/ particular, the porous cylinders transit into several hollow vesicles at φ P = 0.35, the mixtures of hollow vesicle and cylinders transit at φ P = 0.40, 0.45 and distorted cylinders transit at φ P = 0.50. Under weak shear flows, the shear force is insufficiently strong, but porous cylinders can still be broken up by the shear flow. The broken cylinders transit into the hollow vesicles at low concentrations. However, as the phospholipid concentration increases, the number of water molecules participating in the shear flow decreases. This condition leads to the weakening of the shear forces in the cylinders, whilst some cylinders maintain their original shapes, as shown in Fig. 4.
As shear flow increases to moderate strength, such as  γ τ − / 1 = 0.04, the phospholipid cylinders transit into the broken lamellae mixed with the vesicles at the low φ P values of 0.35 and 0.40, as shown in Fig. 5. At the high φ P values of 0.45 and 0.50, the vesicles disappear, and the mixed structure is converted into irregular lamellae. This phenomenon is similar to the shear effects on the other system in which the micelles are elongated along the shear direction 42,58 . These irregular or broken lamellae have bilayered structures due to the amphiphilicity of the phospholipid molecules in the aqueous solutions 59 . The previous experiment has developed an array device to form the multilayer lipid cylinders under shear stress 60 , which differs from the presently reported mixed structures or lamellae but with similar bilayered structures.
For a strong γ τ − / 1  of 0.06, the vesicles develop into lamellae that connect the other parts and then form irregular lamellae in the entire system at the low φ P values of 0.35 and 0.40, as shown in Fig. 6. At the high φ P values of 0.45 and 0.50, a large number of phospholipid molecules participate in the assembly, and the irregular lamellae become planar structures, as shown in Fig. 6. These planar bilayer structures have various spatial separations because the phospholipid concentrations vary in the systems. The bilayer structures possess a large number of flat surfaces when φ P increases. The shear-induced effects observed here are similar to the transition sequences in a lipid monolayer under shear flow by increasing shear rates 61 . However, shear flow does not affect the phospholipid cylinders by transiting into other structures but only sharpens their shapes along the pore directions regardless of shear rates (Fig. S2 in the Supplementary Materials). It is similar to those in applying a thin layer of plane shear, where the distance between the layers is more compact 23 .
Dynamic processes. In this subsection, we focus on the dynamic process of the phospholipid microstructures. We also investigate the average radius of gyration and velocity of the phospholipid molecules with various concentrations and shear rates in these dynamic processes.
We initially concentrate on the dynamic processes of the formation of the hollow vesicle, broken lamella-vesicle mixture and irregular lamella at  γ τ − / 1 = 0.02, 0.04, 0.06 for φ P = 0.35, respectively. Then, we plot the system energies as functions of time step, as shown in Fig. 7. In the shear flow of / 1 γ τ −  = 0.02 (Fig. 7a), the system spends approximately 2200 τ to equilibrate the state, and three stages are observed. The data show that 〈E Tot /k B T〉 = 7.602 www.nature.com/scientificreports www.nature.com/scientificreports/ exists in stable stages. In the first stage, the porous cylinder shape is distorted and pulled by shear forces, and the interactions between the phospholipid and water molecules play a basic role in gathering lipid molecules. The shape distortion stage is maintained at approximately 1050 τ, and the system evolves into a shape adjustment stage, in which the distorted PC-I is deformed. The adjustment stage is maintained approximately 1150 τ, as shown in Fig. 7a. The third stage is the vesicle formation stage, in which vesicles are formed, and the system enters an equilibrium state with stable energy. The dynamic processes in the system are similar to those observed in copolymer systems 13,58,62 . The previous experiments and simulations have reported the dynamic path of spherical vesicles for triblock copolymers in solutions by adding solvent rates 63 . Here, we report a distinct dynamic path from porous cylinders to hollow vesicles for phospholipid polymers in solutions by applying shear flows. Figure 7b,c show the formation processes of the broken lamella-vesicle mixture and the irregular lamella at / 1  γ τ − = 0.04, 0.06 for φ P = 0.35, respectively. The dynamic processes are similar to those observed at / 1  γ τ − = 0.02 for φ P = 0.35, in which the porous cylinder undergoes three stages in the entire process. However, the distortion stage is shortened as the shear rate is increased, and the system spends less time to converge into the formation stages with stable energy, as shown in Fig. 7b,c. This finding indicates that the shear effects on the dynamic processes can be achieved by shortening the convergence time.
We know that the radius of gyration tensor can more accurately characterize the structure and the square radius of gyration tensor R g 2 which is defined as follows 64 : { , , } α β ∈ , N and r i,x are the number of chains and x-coordinate of ith bead, respectively, whilst r c is the mass centre. Based on the radius of gyration tensor, we first discussed the average radius of gyration 〈R g 〉 for the structure, which is defined as follow: so that we investigate the average radius of gyration and velocity for the phospholipid molecules with various concentrations and shear rates in dynamic processes, as shown in Figs 8 and 9. We plot the average radius of gyration 〈R g 〉 for φ P = 0.35 at various shear flows, as shown in Fig. 8a. The shear rate  γ τ − / 1 = 0.0 is also shown in Fig. S3 in the Supplementary Materials. In general, 〈R g 〉 strongly fluctuates in the dynamic processes at γ τ − / 1  = 0.02, 0.04, 0.06. This condition is probably due to the instability of the microstructure in the shear flows because of the movement of cylinders driven by the shear forces; this phenomenon is the same as the behaviour observed in the phospholipid vesicle under shear flows 28,32 . Indeed, the porous cylinder structures transit into the hollow vesicle at γ τ − / 1 = 0.02, in which the shear flows strongly affect their shapes. Porous cylinder transits into the broken lamella-vesicle mixture at γ τ − / 1 = 0.04 and into irregular lamella at / 1 γ τ −  = 0.06, in which the lamellar structure has relative stability under shear flows. As a result, fluctuations cannot be easily observed in the strong shear flows. However, the differences between the average radius of gyration 〈R g 〉 over the entire process are evident due to strong shear effects. Particularly, 〈R g 〉 = 2.94 r c , 3.19 r c , 3.38 r c , as obtained by the linear fittings, which indicate the occurrence of shear effects under strong shear flows. In demonstrating the shear effect clearly, we plot the average radius of the gyration diagonal components in the dynamic processes for φ P = 0.35 and / 1 γ τ −  = 0.06, as shown in Fig. 8b. Three components, namely, 〈R gxx 〉, 〈R gyy 〉 and 〈R gzz 〉, fluctuate strongly; however, 〈R gzz 〉 is larger than 〈R gxx 〉 and 〈R gyy 〉, and 〈R gxx 〉 is nearly equal to 〈R gyy 〉. This result indicates nearly isotropic conformation distributions in the x-y plane in the irregular lamellae. However, the shear force elongates the chains along the z-direction. As a result, 〈R gzz 〉 is considerably larger than 〈R gxx 〉 and 〈R gyy 〉. Therefore, we plot 〈ν x 〉 as a function of time at φ P = 0.35 under / 1 γ τ −  = 0.02, 0.04, 0.06, as shown in Fig. 8c. 〈ν x 〉 increases with time firstly and then elevates into the saturated values of 0.244, 0.605 and 0.973 at  γ τ − / 1 = 0.02, 0.04, 0.06. Evidently, a large shear flow results in a large saturated velocity.
In comparing the effects of phospholipid molecular concentration and shear rate on self-assembled structures, we also plot the average radius of gyration 〈R g 〉 for φ P = 0.50 at various shear flows, as shown in Fig. 9a. The shear rate / 1  γ τ − = 0.0 is also shown in Fig. S4 in the Supplementary Materials. In general, 〈R g 〉 strongly fluctuates in the dynamic processes at γ τ − / 1  = 0.02, 0.04, 0.06, which is consistent with our previous analysis. 〈R g 〉 = 2.95 r c , 3.25 r c , 3.55 r c are also obtained by the linear fittings, which obviously exhibit shear effects under strong shear flows. Compared with the concentration of 0.35, the average radius of gyration 〈R g 〉 has significantly increased because www.nature.com/scientificreports www.nature.com/scientificreports/ the concentration of the phospholipids is higher, more phospholipid molecules participate in the self-assembly and the average radius of gyration increases with concentration. We also studied the three components of the average radius of gyration in the dynamic process with γ τ − / 1 = 0.06, as shown in Fig. 9b. 〈R gzz 〉 is larger than 〈R gxx 〉 and 〈R gyy 〉 , and 〈R gxx 〉 is nearly equal to 〈R gyy 〉. Thus, we can know from the above results that strong shear flows can induce the crystallisation of phospholipid polymers; this phenomenon is consistent with what we observed for other polymers in shear flows in the literature 62 . Finally, we plot 〈ν x 〉 as a function of time at φ P = 0.50 under  γ τ − / 1 = 0.02, 0.04, 0.06, as shown in Fig. 9c. The trend of the velocity distribution at the concentration of 0.5 is the same as the trend at the concentration of 0.35. 〈ν x 〉 increases with time firstly and then elevates into the saturated values of 0.155, 0.27 and 0.74 at γ τ − / 1  = 0.02, 0.04, 0.06. However, the value at 0.5 concentration is significantly less than the value at 0.35, considering that shear flow is applied into the water molecules, and a high concentration leads to a small number of water molecules, thus resulting in low saturated velocities.
The self-assembled structural changes in the applied flow field can be derived from Figs 8 and 9. In comprehensively describing the information of the polymer chains constituting the structure, we calculate the shape factor derived from the square radius of gyration tensor R g 2 , which can be expressed as follows 65,66 : 2 . We plot the shape factor 〈δ〉 as functions of head bead numbers N HB = 3, N TB = 4 under various shear rates at the phospholipid concentration of φ P = 0.35, as shown in Fig. 10.
The graph in Fig. 10 shows the relationship of the shape factor with simulation time at different shear rates. Shape factor 〈δ〉 decreases firstly and then increases as simulation time increases. According to the previous literature, when the shape factor 〈δ〉 is equal to 1 in 3D space, the structure tends to have a rod-like structure. When the shape factor 〈δ〉 is equal to 0, the polymer chain appears as a spherical structure; when the shape factor 〈δ〉 is equal to 0.5, it appears as a circular structure 67 . As shown in Fig. 10a, at the lowest point of shape factor 〈δ〉, the polymer chain of the system forms a state of cluster, which is close to the change of the ellipsoidal structure. This finding can be attributed to the tail end of the polymer chain, which has a certain flexibility. With the extension of simulation time, the polymer chain slowly expands under the action of shearing and finally tends to reach the stable state. The shape is changed to a circular shape given 〈δ〉 = 0.465. Moreover, in the case of weak shear at γ τ − / 1  = 0.02, a hollow vesicle structure is finally formed. Figure 10b shows that the same polymer chain remains to be agglomerated together in the case of moderate shear / 1 γ τ −  = 0.04. The value of the shape factor 〈δ〉 increases www.nature.com/scientificreports www.nature.com/scientificreports/ gently with the increase of simulation time and finally becomes steady at the shape factor 〈δ〉 = 0.52. This finding indicates that in the medium shear state, in addition to the easy formation of a circular structure, the increase of shear force allows the structure of the polymer chain to expand more. Furthermore, the shape is likely an ellipse, which can be explained by / 1  γ τ − = 0.04, in which the structure is a mixed structure of hollow vesicles and broken www.nature.com/scientificreports www.nature.com/scientificreports/ lamellas. As shown Fig. 10c, the trend of the shape factor also decreases firstly and then increases and tends to be gentle. However, in the case of  γ τ − / 1 = 0.06, the shape factor 〈δ〉 = 0.54 indicates that the polymer is under a strong shear rate. The structure of the chain tends to become more elliptical, which facilitates the formation of an irregular lamella.
By studying the shape factor 〈δ〉 under three shear strengths for the concentration of φ P = 0.35, we can clearly see that as the shear strength increases, the shape factor 〈δ〉 gradually becomes larger, and the structure of the polymer chain is changed from the previous ellipse. Moreover, as the shear rate increases, the shape factor becomes gradually larger at the lowest point. This finding can be attributed to the shear factor that causes the shape factor to operate similar to a circular transition. The spherical shape develops towards the ellipse, and the increase of the shear rate also leads to the change in the rate of the polymer structure. In the strong shear state, the polymer chain can rapidly change from an ellipsoid state to an elliptical state. This finding indicates that the applied shear rate can accelerate the change in shape of the polymer chain. We investigated the shear-induced assembly of phospholipid microstructures, by using DPD simulation based on the CG model in water solutions. We selected the shear-free phospholipid cylinders with various cross-sectional shapes and subsequently observed their shear-induced microstructures and the corresponding dynamic processes at various concentrations by changing γ τ − / 1  . We observed the vesicles, lamellae and their mixtures in the water solutions at various φ P and shear rates. The results show that phospholipid molecules form a lamellar structure under strong shear flows, whereas a mixture of vesicle and lamellae can be observed under moderate shear flows. However, weak shear flows easily form phospholipid vesicles in the solutions. Even if the shear rate is enhanced, only the phospholipid porous cylinders can be observed to become a lamellar structure more quickly, and it is necessary to observe the new structure by changing the model. Extending the simulation time does not change the results. www.nature.com/scientificreports www.nature.com/scientificreports/ The dynamic processes indicate that porous cylinders undergo three stages, namely, distortion, adjustment and formation, in the entire dynamic process. Strong shear flows can shorten the distortion and adjustment stages, thereby resulting in the rapid progress into the formation stage. The average radius of gyration and velocity in the dynamic processes indicates that shear rate plays an important role in pulling phospholipid molecules into the solutions under various concentrations. In obtaining more information about the polymer, we also studied the shape factor 〈δ〉 of the polymer chain. The results show that the shape factor 〈δ〉 increases with the increase of simulation time and finally becomes flat. Under different shear rates, the shape factor 〈δ〉 also increases with the increase of shear strength. According to the trend of the shape factor 〈δ〉 under the action of strong shears, the polymer structure can quickly change its initial structure to an ellipsoidal structure, and when the shear force increases, the shape factor 〈δ〉 of the polymer chain is also increased. This work offers insights into the shear-induced microstructures and the dynamic processes of phospholipid cylinders.