Massive quantum regions for simulations on bio-nanomaterials: synthetic ferritin nanocages

Abstract

QM/MM molecular dynamics simulations on the 4His-ΔC* protein cage have been performed using multiple active zones (up to 86 quantum regions). The regulation and nanocage stability exerted by the divalent transition metal ions in the monomer-to-cage conversion have been understood by comparing high level quantum trajectories obtained using Cu²⁺ and Ni²⁺ coordination ions.

---

Biomacromolecules with highly ordered architectures, which play an important role in living organisms, have inspired the creation of new biomaterials with diverse functionalities. Among them, synthetic protein cages (PCs) have received increasing attention since the multi-functionality that each cage can hold is compatible with a wide range of applications.¹ PCs can transport important drugs in the interior voids, while the exterior surface can protect the interior load or can be modified with functional ligands for targeting goals. These distinctive features have made modified PCs potentially interesting vehicles for drug delivery.²

Ferritin PCs are particularly attractive because their subunits self-assemble into spherical nanoparticles with outer and inner diameters of 12 and 8 nm, respectively.³ Ferritin nanocages, which remain intact under non-native conditions (e.g. high pH and temperature), are broken down into individual subunits under acidic conditions, recomposing at neutral pH.³ Due to this phenomenon, ferritin has been proposed as a carrier to transport encapsulated drugs through the human body.⁴

Modified ferritin PCs are highly cooperative symmetrical structures made from a small number of protein building blocks.⁵ Huard et al.⁶ recently engineered a PC, named 4His-ΔC* (Fig. 1a), by self-assembling a recombinant ferritin variant. Although 4His-ΔC* is indistinguishable from the wild-type protein in terms of assembly, stability and function, its monomer-to-cage conversion is selectively regulated by Cu²⁺ at pH 7 (i.e. other divalent ions, such as Zn²⁺ and Ni²⁺, produce a much lower yield of 24-mer formation). Accordingly, 4His-ΔC* opens the possibility of hiding therapeutic agents inside modified cages under physiological conditions for targeting delivery.^7,8 However, the roles played by metal–ligands in the assembly and disassembly of cages, which are relevant for future developments, remain unclear. Explanatory simulations of this complex system require accurate description of the metal ligand associated with each subunit and consideration of the synergies among the 24-subunits. This can be successfully achieved using a multiscale strategy based on the maz-QM/MM-MD scheme, whose QM region is defined as the sum of several disjoint QM subregions (or active zones, AZs) but following the general QM/MM approach to obtain the energy and its gradients (forces). Thus, any particle within an AZ is subject to QM forces from the electronic wave function. However, the interactions with the other AZs are treated through electrostatic interactions (forces from point charges), similar to the rest of the MM region (see the ESI†).^9,10

Fig. 1 (a) Top view of the Cu-4His-ΔC* cage. (b) Alternative views of the Cu-4His-ΔC* cage: C₄ pore interface view through the C₂ symmetry axis (left); C₃ pore interface view through the C₃ symmetry axis (middle); and C₄ pore interface view through the C₄ symmetry axis (right). (c) The general view of the quantum regions in the C₂ dimer interface. Details of the two Cu²⁺ coordination sites through His residues (QM_{C2_4H}), the two Ca²⁺ coordination sites (QM_{C2_Ca}), and the two Cu²⁺ coordination sites at the centre of the monomer (QM_cnt). (d) Details of the quantum region in the C₃ pore interface and (e) the C₄ pore interface in the Cu-cage. The C₂, C₃ and C₄ pore interfaces are highlighted in green, pink and yellow, respectively. The transition metal (Cu²⁺ or Ni²⁺) and the Ca²⁺ ions are highlighted in blue and orange, respectively.

In this work we have used the maz-QM/MM-MD method to examine the stability of both the C₂ dimer interface and the complete cage of 4His-ΔC*. The former consists of two four-helix bundle monomers linked by four metallic coordination sites involving residues of both monomers, while the latter is made of 24 four-helix bundle monomers arranged in octahedral (432) symmetry. The influence of the coordination ion on the stability of the dimer interface has been analysed considering Cu²⁺ and Ni²⁺ coordination ions. Each system, hereafter denoted as Cu- and Ni-dimers, involved six coordination environments (quantum regions) with a total of 188 quantum atoms (Fig. 1b and Fig. S2, ESI†): two coordination sites involving Cu²⁺ or Ni²⁺ transition metal ions (QM_{C2_4H}), two sites coordinating Ca²⁺ ions (QM_{C2_Ca}), and two coordinated metal regions in the center of each monomer (QM_cnt). Besides, the study on the stability of the whole 4His-ΔC* sphere was conducted considering 86 quantum regions (i.e., all coordination spheres of either Cu²⁺ or Ni²⁺, and Ca²⁺ ions in the PC) with a total amount of 2830 quantum atoms. The resulting model (Fig. 1c), denoted as the Cu-cage, comprised 12 C₂ dimer interfaces, 8 C₃ pore interfaces, and 6 C₄ pore interfaces. The average time used for each MD step was about 270 s using 256 processors for the dimer and 1900 s using 928 processors for the whole protein. The present work is the most ambitious and mature application of the PUPIL¹¹ software so far. To the best of our knowledge, this is the first attempt to face a simulation with such a massive number of quantum regions within the QM/MM framework and, therefore, results have provided relevant information on both the 4His-ΔC* PC and the viability of the maz-QM/MM-MD approach when applied to huge biological systems.

Analysis of the maz-QM/MM-MD trajectory of the Cu-dimer in explicit water indicates that the distances between the metal ions and the atoms from the coordination groups are in good agreement with X-ray crystallographic values for the Cu–ferritin cage, though the former are more homogeneous than the latter (Table S1, ESI†). However, metal⋯metal distances between closer AZs experience a significant increase with respect to experimental data. For example, Cu²⁺⋯Cu²⁺ from QM_{C2_4H} regions and Ca²⁺⋯Ca²⁺ from QM_{C2_Ca} are overestimated by 0.26 Å (3%) and 1.06 Å (13%) with respect to experimental data. The largest difference corresponds to the Cu²⁺⋯Cu²⁺ distances of atoms located at the center of each monomer in QM_Cnt regions, which are 0.80 and 0.94 Å (23.5% and 27.6%, respectively) larger than the equivalent X-ray values. These differences are due to the highly solvated environment of the simulated system when compared with experimental data derived from hydrated protein crystals. However, we want to note the shorter metal–ligand distances of QM_Cnt AZ obtained from maz-QM/MM-MD calculations considering the solvated environment that accentuates the separation between the two metals.

On the other hand, metal–ligand distances from simulations on the Ni-dimer are on average larger than those obtained for the Cu-dimer in a range of about 0.5–2% depending of the AZ (Table S1, ESI†). Differences are particularly acute in one of the two QM_{C2_Ca} regions in which Ca²⁺ experiences a change in the coordination pattern (i.e. from three protein ligands in the Cu-dimer to two protein ligands in the Ni-dimer). Besides, metal–metal distances in the QM_{C2_4H} and QM_{C2_Ca} regions are larger for the Ni-dimer than for the Cu-dimer by 5.3% and 8.9%, respectively, whereas the AZs of QM_Cnt shrink by 1% and 3.5%. These geometric differences reflect the different metal–ligand behavior between Cu²⁺ and Ni²⁺ and support the hypothesis that Ni²⁺ introduces changes in the stability of the C₂ interface in 4His-ΔC*.

Table 1 describes the coordination of water molecules to metal ions in Cu- and Ni-dimers, which were derived from the radial distribution functions (RDFs) calculated for metal–O(water) atom pairs. Interestingly, no significant difference is detected between Cu²⁺ and Ni²⁺ hydration patterns. The square planar organization of the QM_{C2_4H} AZ presents one coordinated water molecule that remains stable during the simulation at a distance of ∼1.9 Å. This Cu²⁺–water coordination site was also observed in the crystal structure, even though the Cu²⁺–O(water) distance was of 2.96 Å. The QM_Cnt AZ shows a fractional coordination number since one of the three existing water molecules is shared between the two metal ions. The three water molecules located in the latter AZ were already found in the crystallographic structure with metal–O(water) distances of 1.88, 1.99 and 2.26 Å. Finally, although Ca²⁺ ions in QM_{C2_Ca} AZ are not explicitly coordinated, they interact with five water molecules that form the first solvation shell. The average metal–O(water) distance is 2.50 Å, while in the crystal structure the closer water is at 3.96 Å.

Table 1 Water molecules coordinated to the metal ions in Cu- and Ni-dimers. The distance, d_max (Å), corresponds to the maximum of the first peak in the metal–O(water) RDF, while CN is the coordination number

Interface		Cu-dimer		Ni-dimer
		d max	CN	d max	CN
QMC2_4H		1.85	1.0	1.87	1.0
QMC2_Ca		2.50	4.8	2.50	4.5
QMCnt	M1	1.86	1.3	1.84	1.5
	M2	1.84	2.5	1.86	2.5

---

Hydrogen bonds involving the Tyr39 and Asn74 residues of Cu-4His-ΔC* were reported to be away from the C₂ interface, even though such residues are located near the inner cage surface.⁶ This hydrogen bonding network, which is involved in the stability of the cage structure and the C₂ interface, is sensitive to the divalent metal. In fact, monomer-to-cage conversion was observed for Cu-4His-ΔC* in the presence of divalent ions (e.g. Cu²⁺, Zn²⁺ and Ni²⁺), whereas the ferritin variant MIC1 (4His-ΔC* with Y39E, N74E and P88A mutations) exhibited a drastic dependence on the metal ion for the formation of cages (i.e. the cage was not formed in the presence of Ni²⁺ while the highest extent of cage formation was observed with Cu²⁺). In the MIC1, hydrogen bonds involving Tyr39 and Asn74 were suppressed from the C₂ interface while Glu39, Glu74 and Asp42 added some additional electrostatic interactions.⁶

Simulations on Cu- and Ni-dimers of the 4His-ΔC* protein have clarified the influence of the metal ion on the interface, which can be difficult to model using a simple classical MD approach due to the complex metal–ligand bond coordination. In addition, maz-M/MM-MD trajectories take into account the influence of environmental effects on the stability of the C₂ interface.

Table S2 (ESI†) sorts by occupation the main hydrogen bonds in both Cu- and Ni-dimers for values over 10% (i.e. percentage of the stored snapshots in which the hydrogen bond is formed). Although the Ni²⁺ instability on the ferritin variant MIC1 seems to occur due to the Y39E and N74E mutations, analysis of hydrogen bonds reflects differences between the metal ions before mutation from 4His-ΔC* to MIC1. Thus, the importance of the Tyr39⋯Asn74 hydrogen bond in the C2 interface stability is much higher for the Ni-dimer (83% and 77% for Y39⋯N74′ and Y39′⋯N74, respectively) than the Cu-dimer (38% and 25% for Y39⋯N74′ and Y39′⋯N74, respectively). This is consistent with the fact that, after mutation, the new Glu39 and Glu74 residues induce major repulsive charge interactions within the C₂ interface of the Ni-dimer compared to that of the Cu-dimer. Besides, Asn74⋯Asp42, Tyr39⋯His67 and Tyr39⋯Lys71 hydrogen bonds contribute to the instability of the cage in the presence of Ni.²

The whole Cu-cage metalloprotein was modelled to determine the main interactions on the rest of the ferritin interfaces: C₄ and C₃ (Fig. 1). The former is responsible for the union of four monomers by means of a Cu²⁺ ion, and the latter for the union of three monomers by means of two Ca²⁺ ions. This simulation, which represented a very challenging task due to the huge amount of AZs involved in the system, should also be considered as a proof of principle of the maz-QM/MM MD methodology for large and complex biological systems. Thus, the 86 AZs include the Cu²⁺ and Ca²⁺ ions contained in all the C₂, C₃ and C₄ interfaces of the cage.

Simulation of the Cu-cage involved a huge computational effort since 63% of AZs (i.e. the 54 zones with coordinated Cu²⁺ ions) have a difficult convergence. Furthermore, the communications between different parts of the program were under stress due not only to the vast size of the simulated system but also due to the differences between AZs (i.e. the computational time varies with the number of atoms and metallic ions). As the algorithm responsible for sharing computational resources for minimizing the execution time plays an important role, the software^11,12 was re-engineered to parallelize the construction and simulation mechanism of each AZ independently using the distributed client–server architecture. A client (QMWorker) was created to build, simulate, and manage the information of each AZ, which is exchanged and updated at each iteration cycle with the main program. QMWorkers were distributed equitably between the different calculation nodes of the system in order to avoid memory problems and accelerate the process. On the other hand, the best performance was achieved when the number of calculation nodes is equal to the number of AZs (QMWorker) for an equivalent load distribution among them.

Table 2 lists the metal–ligand distances and the water coordination number on the C₃ and C₄ pore interfaces, derived from the RDFs calculated for the Cu-cage trajectory. Although the square planar organization of the C₄ pore interface is very similar to the QM_{C2_4H} AZs of the C₂ interface, two coordinated water molecules were observed in the former while the latter presents a single water molecule only due to steric effects. Moreover, the coordination of Cu²⁺ at the C₄ interface is similar to that observed by crystallography (Fig. 2c). However, Ca²⁺ ions at the C₃ interface present some difference in water coordination: one is coordinated by 3 Glu residues and 4 water molecules meanwhile the other is surrounded by 3 Asp residues and 5 water molecules (Fig. 2a). The ligands are arranged in an alternated order with maximum probability distances of 2.37–2.38 Å to the closest carboxylic oxygens of Glu and Asp residues. Indeed, the C₃ interface exhibits three additional water molecules compared with X-ray crystallography (Fig. 2b). These are located among Ca²⁺ ions acting as a bridge between the metal ions and the more distant carboxylic oxygens of the ligands, which increases the stability of the C₃ pore interface.

Table 2 Ligands and water molecules coordinated to the metal ions of the Cu-cage. The distance, d_max (Å), corresponds to the maximum of the first peak in the metal–O(water) RDF, while CN is the coordination number

Interface		M–X	d max	CN
QMC4_4H		Cu^2+–H@N	2.18	4.0
		Cu^2+–WAT@O	1.83	2.0
QMC3_Ca	M1	Ca^2+–E@O1	2.38	3.0
		Ca^2+–WAT@O	2.40	4.2
	M2	Ca^2+–D@O1	2.37	3.0
		Ca^2+–WAT@O	2.41	5.0
		Ca^2+ (1)–Ca^2+ (2)	4.01	1.0

---

Fig. 2 (a) Averaged structure of the C₃ pore interface. The metal ions and their coordination sphere are shown. E and D refer to Asp and Glu residues, respectively. Details about the metal coordination with protein residues and water molecules for the (b) C₃ and (c) C₄ pore interfaces are also shown. All figures were extracted from the simulation of the whole Cu-cage system.

In summary, multi-scale simulations presented in this work prove that the C₂ interface of 4His-ΔC* is responsible for the stability of the system in the presence of divalent metallic ions. Ni²⁺ ions affect the hydrogen bonding network at the C₂ interface, inducing the loss of interactions promoted by Cu²⁺ ions. This explains that 4His-ΔC* mutants do not form cage structures when Cu²⁺ is replaced by Ni²⁺.⁶ The multiscale relaxation of the entire Cu-cage protein, in which all Cu²⁺ and Ca²⁺ ions as well as their coordination spheres have been represented at ab-initio level, corroborates the great potential of the maz-QM/MM-MD methodology to model the chemical synergies induced by the different subsystems in very large biological systems. Results have evidenced the key role of C₃ and C₄ interfaces in the stability of the whole nanocage.

This work was supported by MINECO (MAT2015-69367-R), PRACE for access to Curie TN and HN (France at GENCI@CEA and Mare Nostrum at Spain), RES-BSC (Mare Nostrum at Spain), and “ICREA Academia” (C. A.).

Conflicts of interest

There are no conflicts to declare.

Notes and references

1. M. Uchida, M. T. Klem, M. Allen, P. Suci, M. Flenniken, E. Gillitzer, Z. Varpness, L. O. Liepold, M. Young and T. Douglas, Adv. Mater., 2007, 19, 1025–1042.
2. E. J. Lee, N. K. Lee and I.-S. Kim, Adv. Drug Delivery Rev., 2016, 106(Part A), 157–171.
3. M. Kim, Y. Rho, K. S. Jin, B. Ahn, S. Jung, H. Kim and M. Ree, Biomacromolecules, 2011, 12, 1629–1640.
4. Z. Zhen, W. Tang, H. Chen, X. Lin, T. Todd, G. Wang, T. Cowger, X. Chen and J. Xie, ACS Nano, 2013, 7, 4830–4837.
5. X. Liu and E. C. Theil, Acc. Chem. Res., 2005, 38, 167–175.
6. D. J. E. Huard, K. M. Kane and F. A. Tezcan, Nat. Chem. Biol., 2013, 9, 169–176.
7. Z. Zhen, W. Tang, C. Guo, H. Chen, X. Lin, G. Liu, B. Fei, X. Chen, B. Xu and J. Xie, ACS Nano, 2013, 7, 6988–6996.
8. S. Kim, G. S. Kim, J. Seo, G. Gowri Rangaswamy, I.-S. So, R.-W. Park, B.-H. Lee and I.-S. Kim, Biomacromolecules, 2016, 17, 12–19.
9. J. Torras, Phys. Chem. Chem. Phys., 2015, 17, 9959–9972.
10. J. Torras, B. P. Roberts, G. M. Seabra and S. B. Trickey, in Adv. Protein Chem. Struct. Biol., ed. K.-C. Tatyana, Academic Press, 2015, vol. 100, pp. 1–31.
11. J. Torras, E. Deumens and S. B. Trickey, J. Comput. – Aided Mater. Des., 2006, 13, 201–212.
12. J. Torras, G. d. M. Seabra, E. Deumens, S. B. Trickey and A. E. Roitberg, J. Comput. Chem., 2008, 29, 1564–1573.

---

Footnote

† Electronic supplementary information (ESI) available: Methodology description of maz-QM/MM-MD, the molecular system, the simulation protocol, the averaged distances of metal–residue coordination, the hydrogen bond occupancy and partition scheme for the modelling of a large system. See DOI: 10.1039/c7cc09512k

---