Brownian dynamic simulation study on self-assembly of ERPCs

▶Purpose:

1. Study morphology and diffusion of building blocks in dilute solution.
2. Verify the dimerization of the shell chain and the core chain.
3. Confirm the formation of ERPC.
In this part, we used a molecular simulation package named GALAMOST. ^[1]

Figure1. The representative architecture of shell polypeptide (A) and core polypeptide (B). There are 29 beads in shell polypeptide and 19 beads in core polypeptide, with blue and yellow beads representing C and S components, gold and silver representing P and A components, respectively,

▶Model and Method:

We use a coarse-grained model to represent (EK)₁₀-TAT-LZ_E and LZ_R-ELP, as illustrated schematically in Fig.1, according to the estimated ratio of diameter of two building blocks. The length of the ELP part is about 15 while the length of the leucine zipper is about 14 and the length of the functional area is about 5. In Fig.1, each segment of C arms and P arms were given 7 positive charges and 7 negative charges respectively.

In Brownian dynamics, the equation of motion of each bead in the system is governed by Langevin equation

where v ⃗ is the velocity of the bead, γ is the friction coefficient. We fix γ = 1.0 and constructor of a Brownian NVT thermostat object for this system, in which we set the box as 35x35x35, and set the temperature is 1.

In order to fit a certain rigid amino acid chain, we apply a bond angle potential to the coarse-grained bead chains,

where k is the potential constant and θ₀ is the equilibrium angle. Here we set θ = 180°C for C and P types.

The non-bonded interaction between the coarse-grained beads are represented by Lennard-Jones (LJ) potential,

in which ε and σ are the well depth and the bead diameter, respectively. Here r is the distance between two beads, and in this work r_cut = 3.0σ is the cutoff distance beyond which the LJ potential is equal to 0. We set the σ = 1 and ε = 1 for all particle types, as well as the α for S-S is 3.0 and for Z-Z is 0.5 to fit the actual situation since the α positively correlated with the interaction between the two types of particles.

The long-range coulomb interaction between the C type and P type are represented by Ewald summation theory handled in the reciprocal sum.^[2]

The self-energy term

The polypeptide chains are constructed by connecting adjacent beads together with FENE bond potential,

where k is attractive force strength, r₀ is equilibrium length, r_m is maximum bond length. Here we sent k=100, r_m= 1.5 for ‘C-S’, ‘P-Z’, and for others k=500, r_m= 0.96.

A time step of δt = 0.005 is used, and the total simulation steps are 5x10⁶. All simulations are performed using with GALAMOST ^[3] while the results are showing in VMD.^[4]

▶Result & Discussion:

At first, we need to ensure our model can exhibit the same conformation with experiment. Since our experiment of assembly may be done by two steps, our self-assembly simulation is also conducted by two separate steps in order to reproduce the actual situation. Totally, simulation have three parts.

Study conformation of building blocks in dilute solution.

Before studying the properties of self-assembly, we first simulate the conformation of LZ_R-ELP and (EK)₁₀-TAT-LZ_E. The bond angle was set to 180° for C type and P type to give segment of leucine zipper rigidity, meanwhile we set lower k (k=100) value and larger r_m(r_m=1.5) value for FENE bond potential to increase the flexibility of connection between different types of beads (Fig.2). The segment of leucine zipper has a stronger rigidity thus they can maintain a certain configuration in solution while the segment of ELP is very flexible and exhibits a configuration of random coil in dilute solution. This result shows good anastomosis with CD spectrum and secondary structure prediction of the sequence design, consequently, increasing our simulation credibility.

Figure2. The model of LZ_RM-ELP(A) and (EK)₁₀-TAT-LZ_E(B).

Verify the dimerization between shell chain segment and core domain unit.

The (EK)₁₀-TAT-LZ_E and LZ_R-ELP were connected by dimerization of leucine zippers, which is mainly achieved by hydrophobic interaction and electrostatic interaction at certain sites. Our leucine zippers have α-helix structure with some charged residues and hydrophobic residues on different directions at one same side of the helix, besides, the leucine zippers of LZ_R-ELP and (EK)₁₀-TAT-LZ_E have different charges so assembly specificity can be obtained. But in the coarse-grained simulation, we can see the model only have one kind of bead in one side of the chain, so this kind of specificity of dimerization is hard to be accurately reproduced by our model because we can’t make beads anisotropic to let one bead exhibit different properties at different directions. Therefore, we simplified the interaction of leucine zippers and explored whether dimerization of LZ_R-ELP and (EK)₁₀-TAT-LZ_E can be observed in simulation by studying this simplified electrostatic interaction and hydrophilic-hydrophobic interaction from a pair of building blocks.

The C and P segment were given 7 positive charges and 7 negative charges respectively to set up coulomb interaction. The rcut values of Lennard-Jones (LJ) potential for S-C, S-P, S-Z, C-Z, Z-P were set to 1.12 in order to reproduce the effects of hydrophilic-hydrophobic interaction. The simulation result is shown in Fig.3

Figure3. Dimerization between (EK)₁₀-TAT-LZ_E and LZ_R-ELP.

We want to achieve the direction-specific dimerization of the LZ_R-ELP and (EK)₁₀-TAT-LZ_E, namely, to ensure that the ELP and (EK)₁₀ segment stay at different ends of dimer rather than at the same side. ELP segment shows to be hydrophobic while (EK)₁₀ is hydrophilic at 37°C and areas with same kind of water affinity tends to be gathered at the same area. So we used this property of ELPs at 37°C to compensate for the loss of directional specificity caused by coarse granulation, and correct dimer can be observed.

Simulate the formation of ERPC.

Figure4. Assembly of LZ_R-ELP.

In our design, assembly region and functional region are separate to (EK)₁₀-TAT-LZ_E and LZ_R-ELP. And TEM result indicate that LZ_R-ELP can be assembled by itself. Thus, the simulation of ERPCs’ assembly is conducted by assembly of LZ_R-ELP itself and assembly between (EK)₁₀-TAT-LZ_E and LZ_R-ELP.

Therefore, we set the α value of Lennard-Jones (LJ) potential for S-S to 2.0, while the α value for C-C is set to 0.5. 80 chains are randomly generated in a box of 35x35x35. Result was shown in Fig.4, which reveals that micelle with relatively uniform particle size is formed, and it is consistent with the results of experiment part.

For assembly between L(EK)₁₀-TAT-LZ_E and LZ_R-ELP, since the dimerization of leucine zipper is regulated mainly by hydrophobic and electrostatic interaction while hydrogen bonding also plays a supporting role, the interaction between leucine zipper can be strong and complex. To make it feasible to simulate LZ_R-ELP as the core domain and (EK)₁₀-TAT-LZ_E as the "shell" with our coarse-grained model, we connect the (EK)₁₀-TAT-LZ_E and LZ_R-ELP covalently to simulate the assembly of ERPCs, based on the simulation result of Fig. 3, which indicated that the dimerization of the leucine zipper can be achieved and the fact that hydrogen bond cannot be accurately reproduced in our coarse-grained Brownian dynamic simulation.

C segment and P segment are connected through 6 FENE bonds, 40 building blocks are randomly generated in the box under the same conditions (Fig.5). The result demonstrates that the building block forms a spherical micelle with (EK)₁₀ at the outermost, ELPs in inner part.

Figure5. The assembly of ERPC.

In sum, we simulated the process of generating (EK)₁₀-TAT-LZ_E/LZ_R-ELP amphiphilic binary complex and the assembly of ERPCs using Brownian dynamic, and the result is highly consistent with the experiment.

▶Reference:

[1] Zhu YL, Liu H, Li ZW, Qian HJ, Milano G, Lu ZY. GALAMOST: GPU-accelerated large-scale molecular simulation toolkit. J. Comput. Chem. 2013, 34 (25): 2197-2211.
[2] Beckers JVL., Lowe CP, De Leeuw SW. An Iterative PPPM Method for Simulating Coulombic Systems on Distributed Memory Parallel Computers, Molecular Simulation. 1998, 20(6): 369-383
[3] GALAMOST web page: http://galamost.ciac.jl.cn/
[4] Humphrey W, Dalke A, Schulten K. VMD - Visual Molecular Dynamics, J. Mol. Graph. 1996, 14(1): 33-38,27-28.