Yali Wanga and
Xuehao He*abc
aDepartment of Chemistry, School of Science, Tianjin University, Tianjin 300350, China. E-mail: xhhe@tju.edu.cn
bDemonstration Centre for Experimental Chemistry & Chemical Engineering Education, Tianjin University, Tianjin 300350, China
cNational Virtual Simulation Experimental Teaching Centre for Chemistry & Chemical Engineering Education, Tianjin University, Tianjin 300350, China
First published on 10th April 2018
To mimic the unique properties of capsid (protein shell of a virus), we performed Brownian dynamics simulations of the self-assembly of amphiphilic truncated cone particles with anisotropic interactions. The particle shape of a truncated cone in our simulations depended on the cone angle θ, truncated height hc and particle type (AxBy and BxAyBz). The hydrophobic A moieties and hydrophilic B moieties are responsible for attractive and repulsive interactions, respectively. By varying the particle shape, truncated cones can assemble into hollow and vesicle-like clusters with a specific cluster size N. To assemble into hollow vesicles, the truncated height hc must be below a critical value. When hc exceeds this critical value, malformation will occur. The dynamics shows that the vesicle formation occurs in three stages: initially the growth is slow, then rapid, and finally it slows down. The truncated height hc has a stronger impact on the growth kinetics than the cone angle θ or the particle type. We explored how the cluster packing depended on the cooling rate and particle number as well as discussing the relationship between the cluster geometry and the interparticle interactions. Further, we also discuss possible methods to experimentally prepare the truncated cones. The results of our work deepen our understanding of the self-assembly behavior of truncated cones and our results will aid the effective design of particle building blocks for novel nanostructures.
In biological systems, peptides and proteins have the ability to assemble into advanced structures.7 Ionic self-complementary peptides can form nanofibers8 and short peptides can self-assemble into semiflexible nanotapes,9 while surfactant-like peptides can form nanotubes, nanovesicles10,11 and supramolecular structures.12,13 Recent experiments have demonstrated that recombinant fusion proteins can form an amphiphilic complex to further self-assemble into globular protein vesicles by fine turning the protein concentration or temperature;14 and both vesicle size and vesicle geometry can be well-controlled. Beside vesicles, solid particles and hierarchical supraparticles can also be formed through rational protein design.15 Peptides and proteins were treated as cone-like molecules in above mentioned studies, which inspired us to understand their assembly mechanism in more depth and helped us to more effectively design novel building blocks for assembly. With regard to non-biological systems, several assembly systems such as amphiphilic block copolymers,16 metal-polymer amphiphiles17 and patchy colloidal dumbbells18,19 have often been regarded as a cone-like molecules and they can assemble into various structures (e.g., spherical micelles, cylindrical micelles and vesicles) owing to the inherent curvature of building blocks. Cone-like particles or molecules are typical building blocks for assembly and they generally have a nonsymmetrical shape with regard to their two ends. Although much experimental progress has been made in the self-assembly of cone-shaped molecules and particles, some issues such as assembly pathways and the influence of building blocks' parameters including shape and interparticle interaction are not yet exhaustively understood.20
Compared with experimental methods, computer simulations are powerful tools for exploring the self-assembly of specific-shaped particles, although it is difficult to simultaneously consider the particle shapes and interparticle interactions.21–24 Because of computational limitations, the direct calculation of nanoparticle self-assembly using an all-atom model is not feasible; instead one must consider coarse-grained model systems. The widely adopted coarse-grained model is the Kern–Frenkel model25 of patchy spherical particles with short-ranged directional attraction. The whole particle is treated as a “bead” and the patch is either attached to the bead's surface or part of it. Models of this class are formulated in terms of spherical particles with symmetric, regular patches because the interparticle attractions are restricted to orientations in which two patches face each other. In Kern–Frenkel model, the anisotropy of interactions strongly depends on the number of patches and their size.25–27 Another general model with anisotropic interactions was proposed by Zhang et al.,28 who used a more complex coarse-grained model in which the whole particle was constructed from a large number of coarse-grained beads. This coarse-grained model greatly reduced the overall number of interactions that must be calculated compared with an all-atom model. The same or complementary beads within the patches interacted with either attractive or repulsive interactions. This coarse-grained model is capable of describing particles and patches with irregular shape and size and it can easily be modified to study the assembly of nonspherical particles.29,30 In this work, the simplified coarse-grained model was used to characterize the shape of the truncated cone, and a truncated cone included two types of coarse-grained beads to separately represent hydrophobicity and hydrophilicity.
Gaining a quantitative description of a particle's shape and interparticle interactions is a critical and challenging task in assembly simulation of nonspherical particles. Recently, we developed an effective simulation method for the assembly of nonspherical particle with the help of tabulated potential technology.30 Utilizing this model, the self-assembly of two kinds of cone-like particles with a series of cone angles was studied, and various spherical structures with a precise cluster size were formed.24 In that system it was found that clusters began to form simultaneously and the formed clusters aggregated uncontrollably. Those formed clusters were not hollow structures, which considerably prohibits their potential applications,13,31 especially for drug delivery which requires dispersed hollow clusters in solutions.13 In this study, a truncated cone with an amphiphilic surface was designed. Meanwhile, purely repulsive interactions were introduced for hydrophilic coarse-grained beads, which imitates the hydrophilic parts of amphiphilic molecules and it can avoid the cluster–cluster aggregation. The aim of this work was to reveal how the shapes of truncated cone particles affect the self-assembly process and to explore the relationship between complex particle shapes and simple cluster geometry. We employed a quantitative Brownian dynamics (BD) simulation framework for nonspherical particles that was proposed in our prior work.30 This BD framework applied the tabulated potential and the interpolation techniques, and it has advantages of being highly efficient at calculating anisotropic interparticle interactions. This force field between two interacting assembly particles has been calculated before the BD simulations begin. This BD simulation framework not only ensures the precise description of the particle shape and anisotropic interparticle interactions, and it also significantly improves the simulation's efficiency.
Herein, we employ the novel BD simulation framework to study the self-assembly of a series of truncated cones with different particle shapes. General descriptions involved in coarse-grained model of the truncated cones and detailed information on the simulation method are introduced in section 2. The influence of the particle shape on assembled structures, cluster size, cluster geometry, aggregation temperature and cluster growth kinetics are discussed in section 3; the dependence of cooling rate and particle number on the self-assembly and the relationship between interparticle interaction and cluster morphology are also discussed in this section. Furthermore, possible synthesis methods for cone-like particles are suggested, and finally a brief summary of this work is presented in the conclusion section.
For hydrophobic A components, the interaction between A beads is Lennard-Jones potential and that is Uij = 4εAA[(σ/rij)12 − (σ/rij)6], while for the hydrophilic B components purely repulsive interactions are used to reflect its hydrophilicity, Uij = εBB(σ/rij)5. The interaction between two unlike A and B beads is the arithmetic mean of like beads, which is
Here, εAA and εBB are interaction strength parameters, σ is interaction range, which is equal to the diameter of a coarse-grained bead (σ = σAA = σBB = σAB = 0.1 nm), and rij is the distance between the centre of two interacting beads. The interaction strength parameters are defined to ensure that the A beads are strongly hydrophobic and the B beads are weakly hydrophilic (εAA = 2.0 KJ mol−1 and εBB = 0.2 KJ mol−1).
According to the interactions used above, the interaction between two truncated cones (I and J) is the sum of the interactions between coarse-grained beads within the truncated cones,
(1) |
Here, the truncated cones are treated as rigid bodies and the interparticle interactions are directionally dependent. The computation cost for UIJ is high in every simulation step. Instead of this, we used an alternative anisotropic potential, given by
(2) |
The strength parameter εe, range parameter σe and shape parameter re were calculated using a numerical matching method.30 These three parameters are orientation-dependent and the orientation between two particles is represented with three angles αI ∈ [0,π], αJ ∈ [0,π] and βIJ ∈ [0,π], as shown in Fig. 1b. Each orientation angle is divided into 99 equal parts when calculating the tabulated potentials.
The three interaction parameters in eqn (2) can be solved only when eqn (1) yields a potential curve with a negative local minimum.30 However, in the present system UIJ maybe positive without local minimum due to the introduction of the purely repulsive interaction between B beads. The potential fitting method must be modified accordingly, and thus we added an additional attractive interaction Uadd = −1.0/(rIJ − 0.1) to UIJ in eqn (1) in the potential calculation. The additional attraction Uadd can lead to the existence of a negative local minimum in new potential curves at all orientation directions. Consequently, the UIJ in eqn (1) can be well fitted with the anisotropic potential in eqn (2). Note that the additional attractive potential Uadd is excluded in the calculation of translation force and torque in BD simulations. Examples of potential calculation are shown in Fig. S1.†
Because there is a short-ranged interparticle interaction expressed in eqn (2), and to further improve the potential calculation efficiency and to increase the calculation convergence rate, a soft-core potential is introduced when rIJ < rm and the potential calculation is paused at rIJ > rm. Thus, the final form of the anisotropic interparticle interactions can be written as
(3) |
In present coarse-grained model, we only consider the global interactions between particles is attractive or repulsive. Such interaction does not correspond to any real interaction type. In our BD simulation, hydrodynamic interactions is also ignored and the random force represents the randomly fluctuating force exerted on a particle by the surrounding fluid. An implicit solvent model was used here. Our tabulated potential method has an advantage in quantitatively describing specific particle shapes and anisotropic interparticle interactions, unlike other methods such as the patch model which presents anisotropy via the patch size and patch position. The direction-dependent interactions between particles arise from all attractive and/or repulsive beads that make up the particles, not from the “atoms” within the patches. This potential computation method is not only suitable for symmetric, regular particle shape but in principle also for arbitrary shapes. In a BD simulation, the computation cost is entirely unrelated to the number of coarse-grained beads within a particle because the interparticle interactions in BD simulations are calculated with interpolation methods based on tabulated potentials.
To study the self-assembly of truncated cones, a rigid body BD simulation in the canonical ensemble was applied. The simulated box was cubic with a box side l; periodic boundary conditions were used. At each BD step, each particle attempts to move and rotate according to the following equations of motion:
(4) |
(5) |
All BD simulations were cooled from a high-temperature, initially disorder conformation of M particles to a target temperature with stable structures. The low target temperature of each simulation was chosen to allow sufficient time for the particles to assemble into stable structures with unnoticeable change of the system's total energy (the annealing temperature range for each simulation is listed in Table S1†). Unless otherwise stated, an annealing rate, 2.0 × 105 step per K, was applied to avoid kinetically trapped structures. This cooling rate was proved to be affordably slow and yet also fast and efficient enough, compared with using 1.0 × 105 step per K or 3.0 × 105 step per K to produce stable clusters. The specific heat capacity Cv was used to monitor changes of the system's energy. Our studies were restricted to the dilute solution regime, i.e., M = 100 particles were dispersed in a cubic simulation box with a box side l = 10 nm, mainly to reduce the computational time.
As the truncated height hc increased, the inner diameter of the hollow structures also increased. However, the increase of hc was unfavourable for forming ordered clusters. When the truncated height hc was very large, order packing was difficult to achieve and malformations were observed. There exists a critical border for each set of truncated cone at the same θ to assemble into vesicle (see Fig. 2). For AxBy particles, the critical value of hc is 0.5 (Fig. 2b) and for BxAyBz particles the critical values of hc is 0.6, except for θ = 30° (Fig. 2d). The disparity of the critical value of hc at θ = 30° was caused by the fluctuation of simulation data. A truncated cone was evolved from a complete cone by cutting several layers from its top tip. It is not a smooth evolution of particle shape, which may lead to the error of the statistical data in a certain particle shape. The source of the malformations of A68B26 particles at θ = 25°and hc = 0.6 was the development of small local errors in the growing micelles. The clusters continue to grow when there is a malformation, and finally a malformed structure is generated and it is similar to a correct vesicle (Fig. 3a). For A93B36 and A116B47 particles at θ = 30°, 35°and hc = 0.6, the source of the malformation is the inability of malformed clusters to correct themselves after the occurrence of a malformation (Fig. 3b and c). While for BxAyBz particles, the origin of malformations of B18A75B36 particles at θ = 30° and hc = 0.6 is interactions of incomplete clusters; two incomplete clusters can connect together into a single malformed structure and the malformed structure can continue to attract free particles (Fig. 3e). For B16A40B26 and B26A71B47 particles at θ = 25°, 35° and hc = 0.7 the interparticle attractive interactions are not sufficient to counterbalance the repulsive volume effect in the assembly when hc is large enough (Fig. 3d and f).
Cluster geometries with their corresponding particle shapes and cluster size are shown in Fig. 4b and c. There are chemical similarities between amphiphilic molecules and the truncated cones, because both amphiphilic molecules and the truncated cones consist of hydrophobic and hydrophilic parts. The hydrophobic parts of truncated cones pack closely to form hollow and spheroid structure, and the formed structures resemble vesicles comprised of amphiphilic molecules. However, the spheroid structure does not agree with the prediction results of the “critical packing parameter” (CPP).33 Truncated cones with 1/3 < CPP ≤ 1/2 should assemble into cylinders. The difference between the CPP expected structures and the actual assembled structures was also been observed by other groups.21,22,32 Tsonchev et al.32 argued that the relationship between CPP values and the aggregation shapes is not unique and the CPP method should merely be reviewed as a guide. Spherical structures are energetically preferred because there is always a higher number density of cones or truncated cones for spherical structures than for cylinders.32 Thus, it is clear that the problem is to find the densest structure that can be formed using truncated cones. Additionally, Chen et al.21 suggested that the disparity may be limited to rigid cone systems, because the CPP method is suitable for the soft amphiphiles.
In this work, the hydrophobic A components are responsible for hydrophobic attractive interactions and the attractive interactions lead to effective self-assembly, and the hydrophilic “heads” (i.e., B components) only interact with repulsive interactions and the repulsive interactions are unfavourable to effective self-assembly. When the truncation height hc is larger than the critical hc illustrated in Fig. 2, the hydrophilic area of each particle will obviously increase; the result is that the repulsive force between particles also increases significantly, which is not favourable for the effective self-assembly of spherical structures because the attractions between particles are not large enough to inhibit malformed self-assembly. Thus, the densest structures ought to be the results of the balance between the hydrophobic attractive interactions and hydrophilic repulsive interactions. Here, entropic effects were not considered because our previous work24 suggested that the contributions of entropic effect on self-assembly can be omitted, because the entropic effect is very small compared with the interparticle interactions. Furthermore, truncated cones are regarded as rigid bodies and all simulations were carried out in dilute systems (particle number density of 0.1 nm−3). The CPP packing rule is actually more suitable for predicting the behaviour of amphiphiles with soft tails,17,34 and for soft amphiphiles by increasing the tail length, the cluster geometries can transform from spherical micelles to cylinders to bilayers. It is obvious that our current model of truncated cones is not suitable for the CPP packing rule because our hydrophobic “tails” are rigid.
In addition, it is instructive to compare the packing of our truncated cones to the virus capsid (protein shell of a virus) assembly. Most studies of virus capsid assembly process usually impose symmetry constraint, e.g., spherical or icosahedral,35,36 to the cluster geometry before assembly. Thus, the problem is limited to an optimal packing for a few subunits constrained to the surface of a spherical or icosahedral shell.37,38 Our simulation have no special approximation or limitation with regard to the assembly of the truncated cones. We found some cases for which the assembled vesicles (cluster size N = 34, N = 42 and N = 70) are very similar to several virus capsid protein, i.e., cowpea chlorotic mottle virus39 capsid with T = 3 Caspar and Klug (CK) structure (N = 32), a NωV virus40 capsid with T = 4 CK structure (N = 42) and a polyoma virus41 capsid with T = 7 CK structure (N = 72).
Cv = ∂E/∂T = 1/kBT(〈E2〉 − 〈E〉2). |
The temperature that corresponds to the first maximum Cv value is defined as aggregation temperature, Ta and the influence of particle shape on Ta is shown in Fig. 5. It can be seen that the Ta of AxBy particles is higher than that of BxAyBz particles at the same θ and hc. At the same θ, Ta decreases with increasing hc; furthermore, the decrease of Ta for the BxAyBz particles occurs at an obviously higher rate than for AxBy particles. For the same hc, Ta increases with increasing θ for both AxBy and BxAyBz particles. All these phenomena can be explained by considering the differences in the particle hydrophobicity and hydrophilicity. The stronger the hydrophobicity, the better the effective self-assembly, and the more effective self-assembly leads to a higher Ta.
Fig. 5 The dependence of aggregation temperature Ta on the particle shape. The Ta is defined as the temperature corresponding to the first specific capacity heat Cv peak. |
Both the success of forming discrete clusters and the growth kinetics depend delicately on the aggregation energies of the free truncated cones. The influence of hc on the self-assembly process was investigated for AxBy and BxAyBz particles at θ = 30°, as shown in Fig. 6 (the date of other particle shapes is shown in Fig. S2–S4†). The sharp changes in energy are reflected in the Cv curves. At hc = 0.2, the decline in energy is stairs-like with obvious plateaus. When hc is 0.4 and 0.5, the system' total energy decreases continuously and the stairs-like downward trend of energy is no longer obvious; accordingly, there is only one distinct sharp peak in the Cv curve. Compared with the effects of particle type and cone angle θ, the truncated height hc has a strong impact on system energy fluctuations. With the increase of hc, the hydrophobicity of particles weakens and the hydrophilicity of particles becomes stronger, that is, the repulsive interactions between particles increase with the increment of hc, which is not favourable to minimizing the system energy, and thus the system rapidly enters the growth stage of second cluster. The typical configurations of the particles at various temperatures are shown in Fig. 6. Initially, particles are randomly distributed, followed by the formation of curved partial structures, and then a well-packed cluster is generated; thereafter, the remaining free particles begin to grow to form the second cluster until its completion. In the whole assembly process, clusters are formed one-by-one, not together. The reason is that when the temperature decreases to a critical value, a stable nucleus starts to create and particles aggregate into the nucleus till a complete vesicle forms. When a vesicle is formed, the concentration of the free particles in system decreases. It becomes difficult to form a new stable nucleus. Only when the temperature is further decreased, the second stable nucleus can be formed. So, clusters are formed one-by-one in such annealing process.
The role of particle type in the assembly kinetics at θ = 30° was also studied (Fig. 6, and the data for other shapes is shown in Fig. S2–S4†). Three distinct periods are visible: (1) initial lag, closely followed by (2) rapid growth, and finally (3) growth slow-down. Such growth kinetics traits resemble classical nucleation and growth mechanism of virus capsids.42 Stage (1) clearly defines the time range in which nucleation occurs, i.e., in which particles diffuse and collide freely, thus producing the critical growth nucleus. Once the growth nucleus has been formed, the system enters a period of fast growth and the number of free particle decreases rapidly until a robust cluster has been formed. For the same hc value, the role of particle type in the assembly kinetics is not obvious. However, increasing hc has a significant impact on the growth kinetics: a distinct lag period for the growth of the second vesicle is observable at hc = 0.2. But this distinct lag stage vanishes when hc is 0.4 or 0.5, and the second clusters (red curves) begin to grow immediately once the generation of the first clusters (black curves) is complete. In addition, it was found that there are some fluctuations in the cluster size, and these fluctuations indicate the merger of clusters. The merger of clusters will creates oversize misassembled structures, and they are usually generated from complete vesicles and partial vesicles. Fortunately, these oversize misassembled structures can be corrected immediately and they are different to the malformation mentioned in Fig. 3. In brief, varying particle shape, especially hc, changes the interparticle interaction significantly and plays an important role in the vesicle growth kinetics.
The influence of the annealing rate on the assembled structure, energy evolution, and cluster growth kinetics were explored at θ = 30° and hc = 0.2. Three cooling rates were used including 1 × 105 step per K, 2 × 105 step per K and 3 × 105 step per K. When 100 particles were in a 10 nm cubic box, only the cooling rates 1 × 105 step per K and 2 × 105 step per K were considered (Fig. 7a and b); this is because an annealing rate 2 × 105 step per K has been proved to be slow enough for the formation of the final two clusters and the energy change related to the formation of each cluster was evident. The slower cooling rate (3 × 105 step per K) was more useful for distinguishing the formation of the first and second clusters. However, using this cooling rate tripled the computation time (Fig. 7c and d). Thus, the moderate cooling rate of 2 × 105 step per K was used for almost all simulations in this work. The larger system, which comprised 200 particles in a 12.6 nm cubic box with an equal particle number density was also tested under these conditions (Fig. 7b and c). Using the larger system did not alter the specific cluster size N, but it did accelerate the growth of the second cluster (red line). This demonstrates that the specific cluster sizes N for each particle shape is robust and independent of system's particle number. There is a sharp increase of cluster size of cluster 1 in Fig. 7d, and it is due to the contact of a partial vesicle structure and the complete vesicle cluster 1. Because the interaction of A and B beads is attractive, the partial vesicle structure (some side surfaces with A beads are exposed) can contact with the B beads of the complete vesicle (where only B beads on the bottoms of particles can be exposed). Finally, all aggregation structures should be vesicle structures with the maximum A–A contacts corresponding to the lowest energy of the system.
Footnote |
† Electronic supplementary information (ESI) available: The details of particle shapes, examples of potential calculation and additional data of cluster kinetics and potential energy surface. See DOI: 10.1039/c8ra01100a |
This journal is © The Royal Society of Chemistry 2018 |