Toxic Effect of Fullerene and Its Derivatives upon the Transmembrane β2-Adrenergic Receptors

Numerous experiments have revealed that fullerene (C60) and its derivatives can bind to proteins and affect their biological functions. In this study, we explored the interaction between fullerine and the β2-adrenergic receptor (β2AR). The MD simulation results show that fullerene binds with the extracellular loop 2 (ECL2) and intracellular loop 2 (ICL2) of β2AR through hydrophobic interactions and π–π stacking interactions. In the C60_in1 trajectory, due to the π–π stacking interactions of fullerene molecules with PHE and PRO residues on ICL2, ICL2 completely flipped towards the fullerene direction and the fullerene moved slowly into the lipid membrane. When five fullerene molecules were placed on the extracellular side, they preferred to stack into a stable fullerene cluster (a deformed tetrahedral aggregate), and had almost no effect on the structure of β2AR. The hydroxyl groups of fullerene derivatives (C60(OH)X, X represents the number of hydroxyl groups, X = 4, 8) can form strong hydrogen bonds with the ECL2, helix6, and helix7 of β2AR. The hydroxyl groups firmly grasp the β2AR receptor like several claws, blocking the binding entry of ligands. The simulation results show that fullerene and fullerene derivatives may have a significant effect on the local structure of β2AR, especially the distortion of helix4, but bring about no great changes within the overall structure. It was found that C60 did not compete with ligands for binding sites, but blocked the ligands’ entry into the pocket channel. All the above observations suggest that fullerene and its derivatives exhibit certain cytotoxicity.


Introduction
Since the discovery of fullerene (C 60 ), its functions have attracted the intensive interests of scientists [1,2]. Because of the excellent physical and chemical properties of fullerene, it has become one of the most promising nanomaterials, and has been widely used in the pharmaceutical industry [3][4][5][6]. In vivo, fullerene molecules show a tendency to aggregate into nanoparticles which can cross the cell membranes of brains due to their lipophilic properties. Thus, fullerene has been used for drug delivery in vivo because of its high stability and hollow cage-like structure [7,8]. For instance, the diameter of fullerene is comparable to the size of active centers of some viruses, such as HIV, so that it has been utilized to block the entrances of active viral centers and prevent the viruses obtaining nutrition from cells [9][10][11]. Recently, it has been found that nanomaterials composed of fullerene derivatives could stabilize immune effector cells to inhibit the release of proinflammatory mediators, making them potential candidates for treatment of several diseases including asthma, arthritis and multiple sclerosis [12][13][14][15][16].
In addition to the broad applications of fullerene and its derivatives in the pharmaceutical industry, their potential toxicity also drawn the attention of researchers [17,18]. The exposure of largemouth bass to nC 60 led to lipid peroxidation in the brain and glutathione depletion in their gills, affecting the signal transmission and normal expression of proteins [19,20]. Further experimental studies reported that fullerene nanoparticles destroyed the permeability and integrity of cell membranes as they crossed them [21][22][23][24]. Cell membranes are composed of various lipids and transmembrane proteins responsible for the functions and activities of cells. When the fullerene nanoparticles cross the cell membranes, they not only affect the structures and functions of the lipid membranes, but also the transmembrane proteins. Previous experimental studies also revealed that fullerenes could affect the activity and function of proteins, including potassium channel proteins [25,26], HIV proteases [9][10][11]27,28], serum albumins [29][30][31] and glutathione S-transferases [32][33][34].
Using microsecond molecular dynamics (MD) simulations, Monticelli et al. investigated the interaction of C 70 fullerene with Kv1.2 potassium channel proteins. Their simulated results revealed that the fullerenes acted as blockers to change the conformations and functions of the secondary structures of Kv1.2 [35]. Further MD simulations of fullerenes and their related derivatives indicated that they could bind to hydrophobic sites on the surfaces of proteins through specific interactions such as π-π stacking [36,37]. It is generally understood that such kind of π-π stacking interactions exist between the conjugated surface of C 60 and the hydrophobic residues of PHE, TRP, and TYR with aromatic rings.
Fullerene and its derivatives worked as competitors to the ligands of target proteins. Kraszewski et al. found that extracellular fullerene molecules blocked the entry of potassium channel proteins and prevented potassium ions from entering the channel. The intracellular fullerene molecules could enter the channel from the side of the cytoplasm and become stuck in the cavity formed by a few helices of channel protein [25,35]. A series of experimental studies about the cytotoxicity of fullerene derivative C 60 (OH) n have also been reported [38][39][40]. Although there has been a great deal of evidence for the toxicity of fullerene and its derivatives, the atomistic interaction details of fullerenes and related proteins remain unclear so far.
In this study, we performed MD simulation to investigate the toxic effect of fullerenes and their hydroxyl derivatives on the β 2 -adrenergic receptor (β 2 AR), which serves as an important target for cardiac and asthma drugs, and is an extensively studied model system within the G-protein coupled receptor (GPCR). β 2 AR is the protein that mediates muscle relaxation in bronchi, as well as vasoconstriction and vasodilation, the structure of which is shown in Figure 1. It has been found that fullerene derivatives are able to suppress disease onset and to reverse established disease in murine asthma models by decreasing airway inflammation and bronchoconstriction. Therefore, we speculate that fullerene molecules can interact with β 2 AR for the treatment of asthma with associated inflammation. Moreover, Yamawaki et al [41] found that the derivative C 60 (OH) 24 caused cytotoxic injury or cell death in vascular endothelial cells, indicating that exposure to fullerene could represent a risk for atherosclerosis and ischemic heart disease. Furthermore, Wang et al [42]. found that the ECL2 (extracellular loop 2) and the transmembrane helix 7 (TM7) in β 2 AR formed a hydrophobic bridge upon the entry site of the ligand. Moreover, the extracellular top of the binding pocket was shown to be composed of rich hydrophobic residues such as PHE, TYR and LYS. Inspired by the structural features of β 2 AR, we postulate that C 60 molecules could interact with β 2 AR through hydrophobic interaction. The ligands of β 2 AR are mainly bound by polar hydrogen bonds, so the polar fullerene derivatives may also interact with the β 2 AR. Therefore, it is necessary to study the interaction between β 2 AR and fullerene or its hydroxyl derivatives.
Due to the complexity of the biosystem, it is not practical to investigate experimentally the effects of fullerenes and fullerene derivatives upon β 2 AR. Computer simulation can be a good alternative. Although there have been large quantities of results describing the interactions between C 60 and proteins observed by experimental and computational methods, no research has been reported focusing on the interaction between C 60 and β 2 AR. In this work, we used MD simulation to investigate the interactions of fullerene and its hydroxyl derivatives with β 2 AR, to elucidate their toxic effect on the function of β 2 AR based on the simulated trajectories. We mainly aimed to address the following issues: (1) identification of the binding sites of C 60 , C 60 (OH) 4 , and C 60 (OH) 8 in β 2 AR; (2) the effects of mono-C 60 , five C 60 s, C 60 (OH) 4 , and C 60 (OH) 8 on the structure of β 2 AR; (3) whether or not C 60 and C 60 (OH) X (X = 4, 8) compete with the binding site of β 2 AR; (4) whether C 60 and C 60 (OH) X (X = 4, 8) block protein entry. Due to the complexity of the biosystem, it is not practical to investigate experimentally the effects of fullerenes and fullerene derivatives upon β2AR. Computer simulation can be a good alternative. Although there have been large quantities of results describing the interactions between C60 and proteins observed by experimental and computational methods, no research has been reported focusing on the interaction between C60 and β2AR. In this work, we used MD simulation to investigate the interactions of fullerene and its hydroxyl derivatives with β2AR, to elucidate their toxic effect on the function of β2AR based on the simulated trajectories. We mainly aimed to address the following issues: 1) identification of the binding sites of C60, C60(OH)4, and C60(OH)8 in β2AR; 2) the effects of mono-C60, five C60s, C60(OH)4, and C60(OH)8 on the structure of β2AR; 3) whether or not C60 and C60(OH)X (X = 4, 8) compete with the binding site of β2AR; 4) whether C60 and C60(OH)X (X = 4, 8) block protein entry.

Construction of Fullerene and Derivative Models
The initial structures of C60 fullerene (with radius 3.48 Å) and its hydroxyl derivatives C60(OH)X (X = 4 or 8) were built as shown in Figure 2, where X means the number of hydroxyl groups. For the derivative C60(OH)4 with four hydroxyls, we constructed two representative models C60(OH)4_1 and C60(OH)4_2, with their hydroxyl groups located at different surface positions of C60. In the model C60(OH)4_1, the four hydroxyls were located at the close positions of the fullerene spherical surface, with two or three carbon atoms between the two adjacent hydroxyls. In the other model C60(OH)4_2, the four hydroxyls were evenly distributed spherical fullerenes. For the C60(OH)8 model, there were four hydroxyls close to each other, but the other four hydroxyls were far from each of these four hydroxyls. The structures of fullerene and derivatives were optimized using the B3LYP/6-31G* method in the Gaussian09 software [43]. Then, frequency analysis was carried out to obtain the vibrational force constants for all the bonds and angles involved in the fullerene and derivatives. The RESP charges of all atoms in the fullerenes and derivatives were calculated by fitting the molecular electrostatic potential [44,45]. The van der Waals interactions of pairwise atoms were described using the classical Lennard-Jones potential. The parameters of the cross section s and well depth ε in the Lennard-Jones potential were

Construction of Fullerene and Derivative Models
The initial structures of C 60 fullerene (with radius 3.48 Å) and its hydroxyl derivatives C 60 (OH) X (X = 4 or 8) were built as shown in Figure 2, where X means the number of hydroxyl groups. For the derivative C 60 (OH) 4 with four hydroxyls, we constructed two representative models C 60 (OH) 4 _1 and C 60 (OH) 4 _2, with their hydroxyl groups located at different surface positions of C 60 . In the model C 60 (OH) 4 _1, the four hydroxyls were located at the close positions of the fullerene spherical surface, with two or three carbon atoms between the two adjacent hydroxyls. In the other model C 60 (OH) 4 _2, the four hydroxyls were evenly distributed spherical fullerenes. For the C 60 (OH) 8 model, there were four hydroxyls close to each other, but the other four hydroxyls were far from each of these four hydroxyls. The structures of fullerene and derivatives were optimized using the B3LYP/6-31G* method in the Gaussian09 software [43]. Then, frequency analysis was carried out to obtain the vibrational force constants for all the bonds and angles involved in the fullerene and derivatives. The RESP charges of all atoms in the fullerenes and derivatives were calculated by fitting the molecular electrostatic potential [44,45]. The van der Waals interactions of pairwise atoms were described using the classical Lennard-Jones potential. The parameters of the cross section s and well depth ε in the Lennard-Jones potential were taken from the Amber03 force field, with values of 0.34 nm and 0.36 kJ/mol respectively.

Simulation Systems
The initial structure of β 2 AR was extracted from the crystal structure (PDB ID: 2RH1) deposited in the online Protein Data Bank [46]. The crystal structure of β 2 AR contains an endogenous ligand (−)-Isoproterenol [47,48], which is a beta adrenoreceptor agonist used for the treatment of bradycardia (slow heart rate), heart block, and asthma. The (−)-Isoproterenol in the system is abbreviated as 'IPT'.

Simulation Systems
The initial structure of β2AR was extracted from the crystal structure (PDB ID: 2RH1) deposited in the online Protein Data Bank [46]. The crystal structure of β2AR contains an endogenous ligand (−)-Isoproterenol [47,48], which is a beta adrenoreceptor agonist used for the treatment of bradycardia (slow heart rate), heart block, and asthma. The (−)-Isoproterenol in the system is abbreviated as 'IPT'.
To explore systematically the interaction of the β2AR protein with the fullerene and derivatives, we constructed fifteen systems composed of β2AR protein complexes, DOPC To explore systematically the interaction of the β 2 AR protein with the fullerene and derivatives, we constructed fifteen systems composed of β 2 AR protein complexes, DOPC lipids and water molecules. Details of the components in these complex systems are shown in Table 1.
In these complex systems, the initial positions of fullerenes and derivatives were randomly placed at the extracellular or intracellular sides of lipid membranes, indicated by "ex" or "in" in the notation. Initially, to reduce the interaction between the fullerene molecule or fullerene derivatives and β 2 AR, they were spatially separated at distances (defined as the distance between the center of fullerene or derivative's mass and the β 2 -adrenergic receptor surface) of 1.4 nm for C 60 _ex1 (C 60 _ex2), 0.7 nm for C 60 _ex3 (C 60 _ex4), 0.9 nm in C 60 _in1 (C 60 _in2), 1.0 nm for C 60 _IPT_ex, and 0.8 nm for C 60 (OH) 8 _IPT_ex. In the simulation of C 60 _ex5, the fullerene was in the middle of the ECL1 of β 2 AR and the neighboring lipid membranes. The centroid distance between the fullerene and the lipid membrane corresponded to the membrane thickness at the position of the lipid membrane surface. The fullerene of C 60 _in3 was put on the top of the pocket using the helix3, helix4, helix6, and helix7. The IPT of the system C 60 _IPT_ex, C 60 (OH) 8 _IPT_ex was manually placed at the top of the binding pocket. The specific positions of C 60 , C 60 (OH) 4 , and C 60 (OH) 8 IPT molecules in the corresponding systems are shown in Figure 2. Table 1. Details of the constructed complex systems composed of fullerene, fullerene derivatives and β 2 AR, with the 122 DOPC lipids and different numbers of water molecules.

Notations Components Number of DOPCs/Waters
Pure_β 2 AR Without fullerene or fullerene derivative 122/12274 C 60 _ex1, C 60 _ex2 One C 60 in the extracellular side 122/12274 C 60 _ex3, C 60 _ex4 One C 60 in the extracellular side 122/12202 C 60 _ex5 One C 60 in the extracellular side 122/12113 C 60 _in1, C 60 _in2 One C 60 in the intracellular side 122/12035 C 60 _in3 One C 60 in the intracellular side 122/11542 C 60 (OH) 4 _ex1 One C 60 with 4 OH in the extracellular side 122/11153 C 60 (OH) 4 _ex2 One C 60 with 4 OH in the extracellular side 122/11375 5C 60 _ex Five C 60 s in the extracellular side 122/13899 C 60 (OH) 8 _ex1 One C 60 with 8 OH in the extracellular side 122/12059 C 60 _IPT_ex One C 60 and one isoproterenol in the extracellular side 122/12308 C 60 (OH) 8 _IPT_ex One C 60 with 8 OH and one isoproterenol molecule in the extracellular side 122/12902

Simulation Details
The MD simulations for all systems were carried out using the Amber03 force field in the Gromacs 4.4.5 software [49]. In the simulation, the bonds involving hydrogen atoms were constrained using the LINCS protocol [50] and the non-bonded interactions were updated every five integration steps. The DOPC bilayers, fullerenes, β 2 Ars, and water molecules were coupled to the Berendsen thermostat [51] at 300 K with a relaxation time of 0.1 ps. The surface tension of bilayer liquids was maintained at a magnitude of 440 bar/nm. The pressure perpendicular to the liquid surface was kept at 1.0 bar using Berendsen pressure coupling with a relaxation time of 1.0 ps and a water compressibility of 4.5 × 10 −5 bar −1 . The electrostatic interaction of pairwise atoms was estimated using the Particle Mesh Ewald method, and van der Waals interaction within a cutoff 1.2 nm was taken into account. All the systems were first energy-minimized, then heated to 300 K and equilibrated for 300 ns in the NPT ensembles. All data in figures were analyzed based on the equilibrated trajectories of 300 ns MD simulations. Figure 3 shows the various binding patterns of monomeric fullerenes to the β 2 AR proteins. The C 60 molecules in these plots display the final positions in snapshots extracted from MD trajectories. It can be clearly seen that the binding sites of fullerenes were mainly located near the ECL2 or ICL2, except for the system C 60 _in1. During the MD simulations, the positions of fullerenes in C 60 _ex1/2/3/4/5 and C 60 _in2/3 showed almost no large changes compared to their initial positions. The reason for this might be the fact that the β 2 AR is a polar receptor, so it prevents the nonpolar C 60 from entering. As shown in the MD simulations, the C 60 molecules in the aqueous phase attached quickly to β 2 ARs within 5 ns. Since the binding pockets of β 2 AR are hydrophilic, the C 60 molecules prefer to bind to the exposed ECL2 instead of entering it. Figure 4 shows the structural details of the specific amino acids at the binding sites of β 2 AR in C 60 _ex1/2, C 60 _ex3/4, C 60 _ex5 and C 60 _in2 interacting with C 60 . It can be seen that the C 60 molecules were embedded in the cavities formed by the hydrophobic moieties of polar and non-polar amino acids such as GLN, MET, and PHE of ECL2. The aromatic rings of the TYR residues in C 60 _ex3/4, PHE in C 60 _ex5, and PHE of ICL2 in C 60 _in2 formed the π-π stacking interaction to stabilize the conformations of C 60 . located near the ECL2 or ICL2, except for the system C60_in1. During the MD simulations, the positions of fullerenes in C60_ex1/2/3/4/5 and C60_in2/3 showed almost no large changes compared to their initial positions. The reason for this might be the fact that the β2AR is a polar receptor, so it prevents the nonpolar C60 from entering. As shown in the MD simulations, the C60 molecules in the aqueous phase attached quickly to β2ARs within 5ns. Since the binding pockets of β2AR are hydrophilic, the C60 molecules prefer to bind to the exposed ECL2 instead of entering it. Figure 4 shows the structural details of the specific amino acids at the binding sites of β2AR in C60_ex1/2, C60_ex3/4, C60_ex5 and C60_in2 interacting with C60. It can be seen that the C60 molecules were embedded in the cavities formed by the hydrophobic moieties of polar and non-polar amino acids such as GLN, MET, and PHE of ECL2. The aromatic rings of the TYR residues in C60_ex3/4, PHE in C60_ex5, and PHE of ICL2 in C60_in2 formed the π-π stacking interaction to stabilize the conformations of C60.   In Figure 3, it can be observed that the C60 in the C60_in1 simulation eventually entered the lipid membrane, which differed from the C60_in2 simulation with the same initial structure. To gain insight into the dynamic process of C60 entry into the membrane, we extracted a series of snapshots from the MD trajectory and marked the positions of C60 at different simulation times using different colors, as shown in Figure 5. The illustration shows that C60 gradually squeezed into the lipid membrane during the 300ns simulation. Snapshots of the C 60 _ex1/2, C 60 _ex3/4, C 60 _ex5 and C 60 _in2 systems show that the key amino acids of the ECL1 and ECL2 interact with the monomeric fullerenes at the extracellular and intracellular sides.

Binding Sites of Monomeric Fullerenes upon β 2 AR
In Figure 3, it can be observed that the C 60 in the C 60 _in1 simulation eventually entered the lipid membrane, which differed from the C 60 _in2 simulation with the same initial structure. To gain insight into the dynamic process of C 60 entry into the membrane, we extracted a series of snapshots from the MD trajectory and marked the positions of C 60 at different simulation times using different colors, as shown in Figure 5. The illustration shows that C 60 gradually squeezed into the lipid membrane during the 300 ns simulation. Figure 4. Snapshots of the C60_ex1/2, C60_ex3/4, C60_ex5 and C60_in2 systems show that the key amino acids of the ECL1 and ECL2 interact with the monomeric fullerenes at the extracellular and intracellular sides.
In Figure 3, it can be observed that the C60 in the C60_in1 simulation eventually entered the lipid membrane, which differed from the C60_in2 simulation with the same initial structure. To gain insight into the dynamic process of C60 entry into the membrane, we extracted a series of snapshots from the MD trajectory and marked the positions of C60 at different simulation times using different colors, as shown in Figure 5. The illustration shows that C60 gradually squeezed into the lipid membrane during the 300ns simulation.  Figure 6 shows the positions of fullerenes in the snapshots extracted from the MD trajectory, with the specific residues surrounded. In the snapshots at 1ns and 30 ns, the C60 molecules showed strong packing interaction with the aromatic rings of PRO and PHE residues. Subsequently, the C60 molecules gradually moved away from the ICL2 region, and embedded deeply to interact with the hydrophobic residues in helix4, such the PHE and VAL residues. It can be seen that the loop in the ICL2 region completely flipped after 50 ns, so that the interaction became weak. Finally, the C60 molecule was bound to the  Figure 6 shows the positions of fullerenes in the snapshots extracted from the MD trajectory, with the specific residues surrounded. In the snapshots at 1 ns and 30 ns, the C 60 molecules showed strong packing interaction with the aromatic rings of PRO and PHE residues. Subsequently, the C 60 molecules gradually moved away from the ICL2 region, and embedded deeply to interact with the hydrophobic residues in helix4, such the PHE and VAL residues. It can be seen that the loop in the ICL2 region completely flipped after 50 ns, so that the interaction became weak. Finally, the C 60 molecule was bound to the hydrophobic PHE, TYR, and VAL residues in the helix4, and its conformation changed remarkably.

Binding Sites of Multiple Fullerenes upon β 2 AR
Using a high concentration of fullerenes in experiments, multiple C 60 molecules should be present to interact with the β 2 AR protein. To explore the interaction mechanism, we randomly placed five C 60 molecules at the extracellular side of β 2 AR at the beginning of the MD simulation. The separation between C 60 molecules was beyond 15 Å, and the distances between the centroids of C 60 and the protein surface was about 7 Å. Figure 7 shows the traces of C 60 molecules as well as the final stable conformations. The traced dots indicate that although the initial structures of C 60 were separated far apart from each other, four of them, i.e., C 60 _2, C 60 _3, C 60 _4 and C 60 _5, finally tended to form a C 60 cluster. In particular, the cyan dots reveal that the C 60 _2 moved nearly 25 Å from its initial position and finally formed a stable cluster with other C 60 molecules. C 60 _1 could not move around and was trapped in a local stable conformation due to the obstacle of ECL2.
The simulation result shown in Figure 7 indicates that the multiple C 60 molecules tended to form aggregates at the surface of β 2 AR because of their strong hydrophobic effect. Thus, the hydrophobic surface of C 60 preferred to interact with the hydrophobic surface of β 2 AR, so that the entropy of the whole protein-protein interaction was most favorable. Furthermore, the simulation revealed that the formed C 60 cluster remained at the entrance of the β 2 AR binding site, meaning that the high concentration of C 60 molecules might block the ligand binding to β 2 AR and interrupt the function of β 2 AR. hydrophobic PHE, TYR, and VAL residues in the helix4, and its conformation changed remarkably. Figure 6. Snapshots of the system C60_in1 at different simulation times show that the key amino acids such as the PRO232, PHE233, LYS234, GLN236 and LEU238 interacted with the monomeric fullerenes at the extracellular and intracellular sides.

Binding Sites of Multiple Fullerenes upon β2AR
Using a high concentration of fullerenes in experiments, multiple C60 molecules should be present to interact with the β2AR protein. To explore the interaction mechanism, we randomly placed five C60 molecules at the extracellular side of β2AR at the beginning of the MD simulation. The separation between C60 molecules was beyond 15 Å, and the distances between the centroids of C60 and the protein surface was about 7 Å. Figure 7 shows the traces of C60 molecules as well as the final stable conformations. The traced dots indicate that although the initial structures of C60 were separated far apart from each other, four of them, i.e., C60_2, C60_3, C60_4 and C60_5, finally tended to form a C60 cluster. In particular, the cyan dots reveal that the C60_2 moved nearly 25 Å from its initial position and finally formed a stable cluster with other C60 molecules. C60_1 could not move around and was trapped in a local stable conformation due to the obstacle of ECL2.   The simulation result shown in Figure 7 indicates that the multiple C60 molecules tended to form aggregates at the surface of β2AR because of their strong hydrophobic effect. Thus, the hydrophobic surface of C60 preferred to interact with the hydrophobic surface of β2AR, so that the entropy of the whole protein-protein interaction was most favorable. Furthermore, the simulation revealed that the formed C60 cluster remained at the entrance of the β2AR binding site, meaning that the high concentration of C60 molecules might block the ligand binding to β2AR and interrupt the function of β2AR.

Binding Sites of Fullerene Derivatives upon β 2 AR
Because C 60 is nonpolar, the major interaction of fullerene with β 2 AR is hydrophobic. This does not hold true for the fullerene derivatives, since the derivatives C 60 (OH) X (X = 4  Figure 8 shows the stable interaction patterns of C 60 (OH) 4 _ex1, C 60 (OH) 4 _ex2, and C 60 (OH) 8 with the β 2 AR proteins. It appears that the binding capacity of fullerene derivatives was stronger than that of C 60 molecules. For example, the hydroxyl groups in C 60 (OH) 4 _ex1 formed hydrogen bonds with the GLU and ASN residues in ECL2. In C 60 (OH) 4 _ex2, the hydroxyl groups formed hydrogen bonds with the ASP, GLU, and THR residues, and in C 60 (OH) 8 _ex, the hydroxyl groups interacted with the HIE, LYS, ASP, and CYX residues. The hydroxyl groups in the fullerene derivatives acted like claws and firmly hooked β 2 AR so that the derivatives blocked the entrance of the β 2 AR binding channels. Although C60(OH)4_ex1 and C60(OH)4_ex2 each contain four hydroxyl groups, their interaction patterns with β2AR were not entirely similar. In the first 40 ns, the four hydroxyl groups in C60(OH)4_ex1 formed hydrogen bonds with β2AR. After 50 ns, the hydrogen bonds between C60(OH)4_ex1 and β2AR were broken and the hydroxyl groups pointed toward the aqueous phase. This phenomenon was not observed in the simulation of C60(OH)4_ex2, in which the four hydroxyl groups were located far away from each other. By analyzing the interaction between C60(OH)4_ex2 and β2AR, we found that the strong hydrogen bonds formed between C60(OH)4_ex2 and the ASP164 and GLU152 residues in ECL2 were retained throughout the 300ns simulation, which maintained the stable conformation of C60(OH)4_ex2. Comparison of the MD simulations of C60(OH)4_ex1 and C60(OH)4_ex2 revealed that the position of hydroxyl groups had a great influence on the binding capacity of fullerene derivatives. When the hydroxyl groups on the surface of C60 stayed close to each other, stable hydrogen bond interaction with β2AR was difficult to maintain, due to the steric effect.
Meanwhile, the number of hydroxyl groups might affect the binding capacity of derivatives to β2AR. In order to explore the dependence of binding capacity on the number of hydroxyl groups, we also carried out MD simulation for the C60(OH)8_ex system, as shown in Figure 8. This showed that the C60(OH)8 molecule was deeply embedded in the β2AR. The hydroxyl groups behaved like several claws firmly grasping the β2AR, blocking the binding entry of ligands. Strong hydrogen bonds formed between the hydroxyl groups and each of the residues ASP164 and GLU152, and were not broken during the 300ns simulation time. The difference between C60(OH)4_ex2 and C60(OH)8_ex indicates that increasing the number of hydroxyl groups in fullerenes could enhance their ability to bind with β2AR.

Stability of β2AR with Binding
To explore the effect of fullerene binding, and that of its derivatives, on the stability Although C 60 (OH) 4 _ex1 and C 60 (OH) 4 _ex2 each contain four hydroxyl groups, their interaction patterns with β 2 AR were not entirely similar. In the first 40 ns, the four hydroxyl groups in C 60 (OH) 4 _ex1 formed hydrogen bonds with β 2 AR. After 50 ns, the hydrogen bonds between C 60 (OH) 4 _ex1 and β 2 AR were broken and the hydroxyl groups pointed toward the aqueous phase. This phenomenon was not observed in the simulation of C 60 (OH) 4 _ex2, in which the four hydroxyl groups were located far away from each other. By analyzing the interaction between C 60 (OH) 4 _ex2 and β 2 AR, we found that the strong hydrogen bonds formed between C 60 (OH) 4 _ex2 and the ASP164 and GLU152 residues in ECL2 were retained throughout the 300 ns simulation, which maintained the stable conformation of C 60 (OH) 4 _ex2. Comparison of the MD simulations of C 60 (OH) 4 _ex1 and C 60 (OH) 4 _ex2 revealed that the position of hydroxyl groups had a great influence on the binding capacity of fullerene derivatives. When the hydroxyl groups on the surface of C 60 stayed close to each other, stable hydrogen bond interaction with β 2 AR was difficult to maintain, due to the steric effect.
Meanwhile, the number of hydroxyl groups might affect the binding capacity of derivatives to β 2 AR. In order to explore the dependence of binding capacity on the number of hydroxyl groups, we also carried out MD simulation for the C 60 (OH) 8 _ex system, as shown in Figure 8. This showed that the C 60 (OH) 8 molecule was deeply embedded in the β 2 AR. The hydroxyl groups behaved like several claws firmly grasping the β 2 AR, blocking the binding entry of ligands. Strong hydrogen bonds formed between the hydroxyl groups and each of the residues ASP164 and GLU152, and were not broken during the 300 ns simulation time. The difference between C 60 (OH) 4 _ex2 and C 60 (OH) 8 _ex indicates that increasing the number of hydroxyl groups in fullerenes could enhance their ability to bind with β 2 AR.

Stability of β 2 AR with Binding
To explore the effect of fullerene binding, and that of its derivatives, on the stability of β 2 AR, we calculated the RMSDs of β 2 AR backbones for the eight systems, as shown in Figure 9. The RMSDs of all curves, except the C 60 _in1 system, remained around 0.2 nm in the simulations, consistent with the results simulated by Romo et al [52]. The RMSD values fluctuated around 0.2 nm, indicating that the cores of β 2 AR were extremely rigid and less susceptible to external binding with C 60 and its derivatives. For the system C 60 _in1, the RMSD reached 0.27 nm at 300 ns, a little higher than that of pure β 2 AR. The large RMSD change exhibited by β 2 AR in C 60 _in1 was caused by the hydrophobic interaction of the ICL2 and helix3 with C 60 . In contrast, the RMSD values of C 60 _in3 were relatively lower than those of the pure β 2 AR system. The reason for hislies in the fact that the fullerene molecule in C 60 _in3 was trapped in the cavity formed by helix 3, helix 4, helix 6 and helix 7, which hindered the cooperative movement of helices in β 2 AR. Previous studies reported that β 2 AR could be activated by the cooperative movement of seven helices, and the ligand binding relied on the collaborative movement of seven helices [53][54][55][56]. RMSD change exhibited by β2AR in C60_in1 was caused by the hydrophobic interaction of the ICL2 and helix3 with C60. In contrast, the RMSD values of C60_in3 were relatively lower than those of the pure β2AR system. The reason for hislies in the fact that the fullerene molecule in C60_in3 was trapped in the cavity formed by helix 3, helix 4, helix 6 and helix 7, which hindered the cooperative movement of helices in β2AR. Previous studies reported that β2AR could be activated by the cooperative movement of seven helices, and the ligand binding relied on the collaborative movement of seven helices [53][54][55][56]. We further analyzed the RMSD change of secondary structures of β2AR, especially for the key helix4. Figure 10 clearly shows the average structural changes of helix4 in each system, with respect to the crystal structure as a reference. It can be observed that the helix4 in the pure β2AR showed no remarkable structural change, while those in C60_ex2, C60_ex3, C60_ex4, C60_in2, C60_in3, C60(OH)8_ex, and 5C60_ex changed little. However, the helices in C60_ex1, C60_ex5, C60_in1, C60(OH)4_ex1, and C60(OH)4_ex2 were severely distorted from the crystal structures. Figure 11 shows the RMSD values of helices 4 calculated for the nine systems in Figure 9. The RMSD change of helix4 for C60(OH)4_ex1/2 was larger than that of C60(OH)8_ex. This could be attributed to the strong hydrogen bonding interaction between C60(OH)8 and β2AR in C60(OH)8_ex. In addition, the C60 cluster formed in 5C60_ex had no substantial influence on the stability of helix 4. This suggests that increasing the number of C60 molecules has minimal influence on the stability of β2AR. We further analyzed the RMSD change of secondary structures of β 2 AR, especially for the key helix4. Figure 10 clearly shows the average structural changes of helix4 in each system, with respect to the crystal structure as a reference. It can be observed that the helix4 in the pure β 2 AR showed no remarkable structural change, while those in C 60 _ex2, C 60 _ex3, C 60 _ex4, C 60 _in2, C 60 _in3, C 60 (OH) 8 _ex, and 5C 60 _ex changed little. However, the helices in C 60 _ex1, C 60 _ex5, C 60 _in1, C 60 (OH) 4 _ex1, and C 60 (OH) 4 _ex2 were severely distorted from the crystal structures. Figure 11 shows the RMSD values of helices 4 calculated for the nine systems in Figure 9. The RMSD change of helix4 for C 60 (OH) 4 _ex1/2 was larger than that of C 60 (OH) 8 _ex. This could be attributed to the strong hydrogen bonding interaction between C 60 (OH) 8 and β 2 AR in C 60 (OH) 8 _ex. In addition, the C 60 cluster formed in 5C 60 _ex had no substantial influence on the stability of helix 4. This suggests that increasing the number of C 60 molecules has minimal influence on the stability of β 2 AR.  To uncover the mechanism of the helix4 structural change, we analyzed the conformation evolution of β2AR, as shown in Figure 12. We found that it was closely associated with the conformational change of the key residue PHE193 in ECL2. The PHE193 showed an obvious flip in C60_ex1, C60_in1, and C60(OH)4_ex1. Compared to the RMSD curves in Figure 11, we found that these three systems had large RMSD values for the fluctuation of helix4. This is because the fullerene molecules outside the membrane underwent hydrophobic and π-π stacking interaction with PHE193, and the fullerene molecules in the membrane mainly affected the structure of helix4 by changing the structure of ICL2.  To uncover the mechanism of the helix4 structural change, we analyzed the conformation evolution of β2AR, as shown in Figure 12. We found that it was closely associated with the conformational change of the key residue PHE193 in ECL2. The PHE193 showed an obvious flip in C60_ex1, C60_in1, and C60(OH)4_ex1. Compared to the RMSD curves in Figure 11, we found that these three systems had large RMSD values for the fluctuation of helix4. This is because the fullerene molecules outside the membrane underwent hydrophobic and π-π stacking interaction with PHE193, and the fullerene molecules in the membrane mainly affected the structure of helix4 by changing the structure of ICL2. To uncover the mechanism of the helix4 structural change, we analyzed the conformation evolution of β 2 AR, as shown in Figure 12. We found that it was closely associated with the conformational change of the key residue PHE193 in ECL2. The PHE193 showed an obvious flip in C 60 _ex1, C 60 _in1, and C 60 (OH) 4 _ex1. Compared to the RMSD curves in Figure 11, we found that these three systems had large RMSD values for the fluctuation of helix4. This is because the fullerene molecules outside the membrane underwent hydrophobic and π-π stacking interaction with PHE193, and the fullerene molecules in the membrane mainly affected the structure of helix4 by changing the structure of ICL2. Therefore, it was postulated that the structural changes of helix4 would be directly related to the structural changes of PHE193. Therefore, it was postulated that the structural changes of helix4 would be directly related to the structural changes of PHE193.

Competitive Mechanism of C60/C60(OH)8 and Isoproterenol
Previous experiments have shown that fullerene molecules and fullerene derivatives may compete with ligands for protein binding sites. In our study, we built two systems, labeled as C60_IPT_ex and C60(OH)8_IPT_ex, to investigate the competitive mechanisms of C60 and C60(OH)8. Here 'IPT' is used as the abbreviated name for isoproterenol. Although, the C60 molecule and isoproterenol were both placed above the binding site, we found that C60 acted as a single fullerene and had no influence on the isoproterenol within the simulation time, as shown in Figure 13. This result can easily be understood by the hydrophobic nature of pure C60. Then we placed C60(OH)8 and isoproterenol together, and found that C60(OH)8 quickly combined with β2AR through hydrogen bonds, and the isoproterenol was relegated to the outside. It is interesting that the benzene ring and C60(OH)8 were packed together. When there were many hydroxyl fullerene derivatives outside the cell membrane, they were likely to compete with ligands for the binding site, thus affecting the signal transduction of β2AR to the following proteins. 3.5. Competitive Mechanism of C 60 /C 60 (OH) 8

and Isoproterenol
Previous experiments have shown that fullerene molecules and fullerene derivatives may compete with ligands for protein binding sites. In our study, we built two systems, labeled as C 60 _IPT_ex and C 60 (OH) 8 _IPT_ex, to investigate the competitive mechanisms of C 60 and C 60 (OH) 8 . Here 'IPT' is used as the abbreviated name for isoproterenol. Although, the C 60 molecule and isoproterenol were both placed above the binding site, we found that C 60 acted as a single fullerene and had no influence on the isoproterenol within the simulation time, as shown in Figure 13. This result can easily be understood by the hydrophobic nature of pure C 60 . Then we placed C 60 (OH) 8 and isoproterenol together, and found that C 60 (OH) 8 quickly combined with β 2 AR through hydrogen bonds, and the isoproterenol was relegated to the outside. It is interesting that the benzene ring and C 60 (OH) 8 were packed together. When there were many hydroxyl fullerene derivatives outside the cell membrane, they were likely to compete with ligands for the binding site, thus affecting the signal transduction of β 2 AR to the following proteins.

Conclusions
In this work, we studied for the first time the interactions between fullerene and fullerene derivatives and β2AR, using a series of all-atom MD simulations. We established fifteen systems by placing fullerene and fullerene derivatives at different positions outside and inside the membrane. Our results show that the unmodified fullerene was unable to enter the ligand binding sites, instead mainly binding with the loop of ECL2 and ICL2 on β2AR through hydrophobic interactions and π-π stacking interactions, especially with residue PHE193 on ECL2, and PHE and PRO residues on ICL2. In addition to having a strong hydrophobic interaction with nonpolar amino acid residues, fullerene molecules also underwent obvious hydrophobic interaction with hydrophobic carbon atoms on polar amino acids, such as GLN, GLU, ASP, LYS, TYR, and THR. In most cases (except C60_in1), fullerene and derivatives did not enter the lipid membrane, but remained around the β2AR. In the C60_in1 trajectory, due to the π-π stacking interactions of the fullerene molecule with PHE and PRO residues on ICL2, ICL2 completely flipped towards the fullerene direction along with the movement of fullerene slowly into the lipid membrane. Five fullerenes tended to stack into a stable fullerene cluster (a deformed tetrahedral aggregate) rather than interacting alone with β2AR. In addition, as the simulation time extended, the position of the fullerene aggregates became almost motionless, just above the receptor binding pocket, which would block the movement of ligands into protein channels and hinder the ligands' binding to proteins. The binding capacity between polyhydroxy derivatives and β2AR was stronger compared with unmodified fullerene molecules.
The hydroxyl groups of fullerene derivatives C60(OH)4 and C60(OH)8 can form strong hydrogen bonds with the ECL2, helix6, and helix 7 of 2AR. The hydroxyl groups firmly grasped the β2AR receptor like several claws, blocking the binding entry of ligands. The position and number of hydroxyl groups on fullerene derivatives also seriously affected the interaction between the two groups. The more dispersed the hydroxyl position and the greater the hydroxyl number, the stronger was the binding activity with the β2AR.

Conclusions
In this work, we studied for the first time the interactions between fullerene and fullerene derivatives and β 2 AR, using a series of all-atom MD simulations. We established fifteen systems by placing fullerene and fullerene derivatives at different positions outside and inside the membrane. Our results show that the unmodified fullerene was unable to enter the ligand binding sites, instead mainly binding with the loop of ECL2 and ICL2 on β 2 AR through hydrophobic interactions and π-π stacking interactions, especially with residue PHE193 on ECL2, and PHE and PRO residues on ICL2. In addition to having a strong hydrophobic interaction with nonpolar amino acid residues, fullerene molecules also underwent obvious hydrophobic interaction with hydrophobic carbon atoms on polar amino acids, such as GLN, GLU, ASP, LYS, TYR, and THR. In most cases (except C 60 _in1), fullerene and derivatives did not enter the lipid membrane, but remained around the β 2 AR. In the C 60 _in1 trajectory, due to the π-π stacking interactions of the fullerene molecule with PHE and PRO residues on ICL2, ICL2 completely flipped towards the fullerene direction along with the movement of fullerene slowly into the lipid membrane. Five fullerenes tended to stack into a stable fullerene cluster (a deformed tetrahedral aggregate) rather than interacting alone with β 2 AR. In addition, as the simulation time extended, the position of the fullerene aggregates became almost motionless, just above the receptor binding pocket, which would block the movement of ligands into protein channels and hinder the ligands' binding to proteins. The binding capacity between polyhydroxy derivatives and β 2 AR was stronger compared with unmodified fullerene molecules.
The hydroxyl groups of fullerene derivatives C 60 (OH) 4 and C 60 (OH) 8 can form strong hydrogen bonds with the ECL2, helix6, and helix 7 of 2AR. The hydroxyl groups firmly grasped the β 2 AR receptor like several claws, blocking the binding entry of ligands. The position and number of hydroxyl groups on fullerene derivatives also seriously affected the interaction between the two groups. The more dispersed the hydroxyl position and the greater the hydroxyl number, the stronger was the binding activity with the β 2 AR.
We further explored the effect of fullerene and its derivatives on β 2 AR structure. The results indicated that fullerenes and their derivatives did not have an obvious impact on the whole structure of the receptor, but could affect the structure of helix4 in some systems, especially in the systems C 60 _ex1, C 60 _ex5, C 60 _in1, C 60 _in3, C 60 (OH) 4 _ex1, and C 60 (OH) 4 _ex2. The helix4 structures in the above systems were seriously distorted, with the upper part twisted to the right, and the lower part distorted to the left. Moreover, the more hydroxyl in the fullerene derivative, the smaller was the effect on the β 2 AR structure.
In summary, fullerenes and fullerene derivatives can interact with β 2 AR by hydrophobic interactions, π-π stacking interactions and hydrogen bonding. This combination may block the entry of ligands into protein binding sites, and also may compete with ligands for binding sites, thereby affecting the normal biological function of β 2 AR. Therefore, in future medical research it is important to consider the effect of fullerene and its derivatives on β 2 AR, to further reduce its toxicity to biological cells. Funding: This work was supported by the Natural Science Foundation of Shandong Province (ZR2020MA075, ZR2020QA051), National Natural Science Foundation of China (12104262). We acknowledge the support of the NYU-ECNU center for Computational Chemistry at NYU Shanghai. We also thank the supercomputer center of ECNU for providing computer time.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.