key: cord-325282-20l9xcmg
authors: Helal, Mohamed A.; Shouman, Shaimaa; Abdelwaly, Ahmad; Elmehrath, Ahmed O.; Essawy, Mohamed; Sayed, Shireen M.; Saleh, Amr H.; El-Badri, Nagwa
title: Molecular basis of the potential interaction of SARS-CoV-2 spike protein to CD147 in COVID-19 associated-lymphopenia
date: 2020-09-16
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2020.1822208
sha: 
doc_id: 325282
cord_uid: 20l9xcmg

Lymphopenia is considered one of the most characteristic clinical features of the coronavirus disease 2019 (COVID-19). SARS-CoV-2 infects host cells via the interaction of its spike protein with the human angiotensin-converting enzyme 2 (hACE2) receptor. Since T lymphocytes display a very low expression level of hACE2, a novel receptor might be involved in the entry of SARS-CoV-2 into T cells. The transmembrane glycoprotein CD147 is highly expressed by activated T lymphocytes, and was recently proposed as a probable route for SARS-CoV-2 invasion. To understand the molecular basis of the potential interaction of SARS-CoV-2 to CD147, we have investigated the binding of the viral spike protein to this receptor in-silico. The results showed that this binding is dominated by electrostatic interactions involving residues Arg403, Asn481, and the backbone of Gly502. The overall binding arrangement shows the CD147 C-terminal domain interacting with the spike external subdomain in the grove between the short antiparallel β strands, β1’ and β2’, and the small helix α1’. This proposed interaction was further confirmed using MD simulation and binding free energy calculation. These data contribute to a better understanding of the mechanism of infection of SARS-CoV-2 to T lymphocytes and could provide valuable insights for the rational design of adjuvant treatment for COVID-19. Communicated by Ramaswamy H. Sarma

In December 2019, an outbreak of a pneumonia-causing virus was first reported in Wuhan, China (Lu, Stratton, et al., 2020; Zhou et al., 2020) , and has been spreading globally ever since (Sohrabi et al., 2020) . The virus was found to phylogenetically belong to the severe acute respiratory syndrome coronavirus (SARS-CoV) (Lu, Zhao, et al., 2020) . The virus was named SARS-CoV-2 (C. S. G. of the International, 2020; Gorbalenya, 2020) , and the disease named coronavirus disease-19 . The entry of the virus into host cells is facilitated by binding of its transmembrane spike (S) protein with angiotensin-converting enzyme 2 (ACE-2) receptor (Hoffmann et al., 2020) . The spike protein is a class I viral fusion transmembrane glycoprotein that plays a dual role in viral entry and fusion with host cell membrane (Bosch et al., 2003; Li, 2016) . Cryo-electron microscopy (cryo-EM) revealed the structure of S protein as trimeric spikes composed of three S1 and S2 heterodimer subunits, protruding out of the viral envelope (Gui et al., 2017; Wrapp et al., 2020) . The S1 subunit facilitates viral binding to cell receptor through the receptor-binding domain (RBD), whereas the S2 subunit encloses a hydrophobic fusion domain, allowing viral fusion (Song et al., 2018) . Cleaving the S1-S2 boundary is primed by host cell proteases such as transmembrane protease serine 2 (TMPRSS2), and is required for viral activation and entry (Hoffmann et al., 2020) . TMPRSS2 and ACE-2 are highly expressed in the epithelial tissue of multiple organs and thus significantly contribute to viral pathogenicity, clinical manifestations, and mortality (Sungnak et al., 2020) .

One of the clinical indicators of the severity of SARS-CoV-2 is low lymphocytic count (lymphopenia) , especially in the CD8 þ , CD4 þ , and CD3 þ T lymphocyte populations (Zeng, 2020; Zheng et al., 2020) . Lymphopenia is a risk factor in viral infections, and was reported to predict the severity of the disease in COVID-19 patients (Gui et al., 2017) .

Although several studies examined the causes of lymphopenia during viral infections (Lalueza et al., 2017 (Lalueza et al., , 2019 , the mechanisms underlying this feature are not fully understood in COVID-19 patients. Proposed mechanisms for the observed lymphopenia include direct invasion of lymphocytes , hyperlactic acidemia causing suppression of lymphocyte proliferation (Fischer et al., 2007) , as well as an inflammatory cytokine storm resulting in lymphocyte apoptosis (Liao et al., 2002) , similar to that reported earlier with MERS-CoV infection (Chu et al., 2016) .

Angiotensin-converting enzyme 2, or ACE-2 receptor, has been shown to provide an entry point to SARS-CoV-2 in different types of human cells Zhou et al., 2020) . However, T lymphocytes were shown to be consistently negative for ACE-2 receptors (Hamming et al., 2004) , suggesting an alternate route for viral entry. CD147 was recently proposed as a probable route for SARS-CoV-2 invasion into Vero E6 host cells (Wang, 2020) . CD147 (also known as Basigin or EMMPRIN) is a transmembrane glycoprotein that was originally identified in 1992 as a T lymphocyte activation-associated antigen (Kasinrerk et al., 1992) . Further research confirmed its expression on activated T lymphocytes and its relatively weaker expression in resting T lymphocytes, and that its expression increases rapidly upon T lymphocyte activation (Koch et al., 1999) . CD147 induces matrix metalloproteinase (MMP) production during physiological and pathological events including inflammatory response, wound healing, and tissue homeostasis (Cruzat et al., 2018; Jin et al., 2019; Kato et al., 2009; Ungern-Sternberg et al., 2018) . Interestingly, CD147 receptor is expressed on red blood cells (RBCs), providing a route for malaria entry (Zhang et al., 2018) . It is worth mentioning that PfRH5 complex on the cell membrane of Plasmodium falciparum, the causative agent of malaria, interacts with CD147 on human RBCs, a binding which was necessary for invasion by all strains when tested in vitro (Crosnier et al., 2011) . Anti-CD147 antibodies were shown to prevent entry of P. falciparum and SARS-CoV-2 into host cells (Bian, 2020; Zhang et al., 2018) .

To understand the mechanism of interaction of the SARS-CoV-2 spike protein with the CD147 receptor, we have performed a four-stage in-silico study. First, the binding site of both proteins had to be accurately predicted, taking into consideration the few number of their available crystal structures. Second, molecular docking was performed using three different protein-protein interaction servers. Then, the best docking pose was analyzed in MD simulation to check its stability. Finally, the binding free energy was estimated, using snapshots of the MD simulation, to get in-depth knowledge of the binding determinants.

The study starts with determination of the most probable binding sites for spike and CD147. This was specially challenging for the latter because only few crystal structures were reported for this membrane glycoprotein. Later, the molecular docking was performed using three different web servers and the best pose was selected using the PDBePISA website. Then, the best complex was subjected to a 100 ns MD simulation in NPT ensemble. Finally, free energy calculation was performed to gain insights into the major contributing forces in binding and the critical residues involved.

The SARS-CoV-2 receptor binding domain (RBD) was obtained from its crystal structure complexed with human angiotensin-converting enzyme 2 (hACE2) (PDB ID: 6LZG), while that of CD147 was isolated from its crystal structure complexed with malaria invasion protein, reticulocyte-binding protein homologue (RH5) attained from the parasite Plasmodium falciparum (PDB ID: 4U0Q). The crystal structures of both targeted proteins, SARS-CoV-2 and CD147, were closely examined for integrity and the active residues taking part in the interaction were identified. The surface of CD147 was examined using the Computed Atlas of Surface Topography of Proteins (CASTp) server. The default probe radius used by CASTp is 1.4 Å. We opted to use a smaller probe of 0.8 Å to get more exhaustive coverage of all the potential binding sites on the CD147 surface.

Prior to the docking simulation, both proteins, spike and CD147, were checked for errors and prepared in the Protein Preparation Wizard, within MOE software. The preparation steps include removal of any bound ligands, water of crystallization, and salts. Also, all the hydrogens were added and the protein termini were capped. Both proteins were minimized individually in MOE to a gradient of 0.001 kcal/mol.Å 2 . First, all heavy atoms were kept fixed, and the hydrogen atoms were allowed to optimize their positions. Second, the entire backbone of the protein was kept fixed and the side chains were minimized using the same parameters. Finally, the whole molecule including the backbone was minimized (Helal et al., 2011; Helal & Avery, 2012) . Three different webdocking servers, HADDOCK v2.2, ZDOCK v3.0.2, and HawkDock 2018, were utilized to computationally predict the protein-protein interactions of the two unbound protein structures. The docking trial on each server was repeated twice using the default settings; in the first run no restraints was given to the server so that it can freely determine the binding region (interface) for each protein. This step was performed to find out to which extent the software would bias toward a specific docking pose. In the second run on each server, specific active residues were defined for the software to be considered as restraints (see the discussion part). These residues were selected based on our understanding of the crystal structures of both SARS-CoV-2-RBD/hACE2 and CD147 RH5 as well as the prediction by CASTp server.

Molecular dynamics (MD) simulation was performed using GROMACS 2019 software package and the all-atom OPLS force field. The protein complex was solvated in a cubic box of SPC water model of 100 Â 100 Â 100 Å size and Periodic boundary conditions were implemented. Ionizable residues were assigned standard ionization states at PH 7. The system was then neutralized by the addition of five sodium ions. The total number of atoms in the system, including the water and ions, was 117,013. Short-range interactions were treated with a 10 Å cutoff for both Lennard-Jones and Coulomb potentials, while the Particle Mesh Ewald (PME) algorithm was used for long-range electrostatics. The MD simulation was performed in four steps. During the first three steps, a force constant of 1000 kJ mol À1 nm À2 was used to restrain all the heavy atoms to preserve the original fold of the proteins. First, the solvated protein system was minimized using a 5000 steps of steepest descent algorithm to resolve any steric clashes or inappropriate geometry. Then, to ensure a reasonable starting structure, the system was equilibrated under constant Number of particles, Volume, and Temperature (NVT) ensemble for 100 ps using Berendsen thermostat (Golo & Sha K %tan, 2002) . The second round of equilibration was performed under constant pressure (NPT ensemble) using the Parrinello-Rahman barostat for an additional 100 ps (Tuble et al., 2004) . Finally, the position restraints were released and the system was simulated in the production run under an NPT ensemble (Nos e-Hoover thermostat and Parrinello-Rahman barostat) for 100 ns using a time step of 2 fs. Simulation results were analyzed using the Visual molecular dynamics (VMD) software, ver. 1.9.3 (Humphrey et al., 1996) . Pymol graphical software ver. 2.3 was used for the generation of the figures of protein-protein interactions.

Binding free energy was calculated using the Molecular Mechanics Poisson Boltzmann Surface Area (MM/PBSA) method with the g_mmpbsa script. This method calculates the binding energy using molecular mechanics potential energy (electrostatic and van der Waals interactions) and free energy of solvation (polar and nonpolar solvation energies). A total number of 100 snapshots, extracted at regular intervals from the whole MD trajectory, were used for calculating each energy term. The average for each term was then calculated. In addition, a per-residue binding free energy decomposition at each snapshot was performed and the average binding energy for each residue was calculated and depicted. The solvent-accessible surface area (SASA) model was used to compute non-polar solvation energy with a surface tension of 0.02267 (kJ/mol 2 ) and a probe radius of 1.4 Å.

The recently reported crystal structure of SARS-Cov-2 spike protein complex with ACE2 (PDB ID: 6LZG) reveals that the virus utilizes the external subdomain of the spike Receptor Binding Domain (RBD) to recognize the human ACE2 receptor . This external subdomain is composed mainly of a flexible loop, S438-Y508, connecting two anti-parallel small b strands. Also, a part of this loop, the segment F486-G502, represents the most solvent accessible part of the spike RBD and makes direct contact with ACE2 in the crystal structure. In addition, in a recent in-silico study, it was proposed that a small circular region of this loop, C480-C488, could be responsible for the interaction with the human GRP78 (Ibrahim et al., 2020) . Collectively, these findings have convinced us to use the solvent-accessible F486-G502 segment as the interaction site of the spike protein in our docking simulation. Regarding the CD147 receptor, the only crystal structure available of this protein, bound to another biological molecule, is for its complex with the PfRH5 protein (PDB ID: 4U0Q) (Wright et al., 2014) . In this structure, the PfRH5 binds to the cleft between the two domains of CD147 assisted by a network of H-bond and salt bridge interactions. Notably, the glycoprotein CD147 is a relatively small protein and few pockets could be noticed on its surface. To get an accurate prediction of the potential binding site on CD147, we have performed a comprehensive search for its surface pockets using the Computed Atlas of Surface Topography of Proteins (CASTp) server (Tian et al., 2018) . This tool enables an exhaustive and quantitative characterization of the protein topographic features, providing a complete list of the available pockets with the amino acid residues lining these pockets. Figure 1 (A) shows the pockets with significant volume (>15 Å 3 ) on CD147 surface. Among them, pockets 1, 3, 4, and 6 are the most solvent accessible and constitute a nearly continuous grove extending at the interface of the two domains of CD147. Furthermore, according to the crystal structure 4U0Q, the PfRH5 protein uses pockets 1 (124.67 Å 3 ) and 6 (72.57 Å 3 ) for interaction with CD147. In addition to the CASTp server, we have also used the CPROT tool (Bonvinlab) which is an algorithm for the prediction of protein-protein interface merging six calculation methods into a consensus result (de Vries & Bonvin, 2011) . Interestingly, the CPROT results were very similar to those of CASTp with most residues of pockets 1 and 6 identified (Figure 1(B) ). Therefore, we have selected the residues of these pockets, 1 and 6, as the active residues in our docking study, namely Ser78, Asp79, Gln81, Trp82, Gly103, Pro104, Pro105, Arg106, Val 107, Lys108, Lys127, Ser128, Glu129, Ser193, and Asp194. (Table S1 , Supporting Information).

Molecular docking was performed using three different protein-protein docking servers, namely HADDOCK v2. interacting, residues as the segment F486-G502 of the spike RBD and the residues of pockets 1 and 6 of the CD147 as illustrated in the previous section. In order to evaluate the docking results of the three software servers, the top ranking pose from each was examined using the PDBePISA online tool, developed by The European Bioinformatics Institute, EMBL-EBI, for the exploration of macromolecular interfaces (Krissinel & Henrick, 2007) . Moreover, a Molecular Mechanics energies combined with Generalized Born and Surface Area (MM/GBSA) calculation, as implemented within the Hawkdock server, was used to estimate the binding affinities (Hou et al., 2011; Zhong et al., 2020) . The results of the examination of the top complex generated by each software are illustrated in Table 1 . Interestingly, all the three servers returned very similar arrangements of the interacting proteins with RMSD of 2.23 for ZDOCK/HADDOCK and 0.24 for HawkDOCK/HADDOCK. Moreover, for each server, performing the docking run without specifying any constraints returned quite a similar docking pose, to the one with constraints, within the top ten out of more than one hundred solutions (data not shown). Based on the RMSD values and visual examination, it was noticed that the poses proposed by HADDOCK and HawkDOCK are very close and almost converging. We decided to use the pose generated by the HADDOCK server as it shows a higher interface surface area of 1,069 Å 2 and the best score of À68.00 in the MM/GBSA calculations (Table 1 ).

In addition, thorough investigations were carried out using PDBePISA tool to determine the key interacting residues in the selected docked complex generated by the HADDOCK server. Only interfacing residues within 4 Å van dar Waals distance cutoff between the SARS-CoV-2 RBD and CD147 were examined by the PDBePISA algorithm (Table S3 , Supporting Information). It can be easily noted that the docked SARS-CoV-2 RBD adopted a similar orientation to that observed in its cocrystal structure with hACE2, which highlights the binding role of the region that extends between the residues C480 and Q505 (Figure 2 ). The spike external subdomain forms a surface to which parts of CD147 domains dock. This surface is dominated by the interaction of the C-terminal domain of CD147 which binds to the spike external subdomain in the grove between the short antiparallel b strands b1' and b2' and the small helix a1' (Figure 2(A) ). This grove makes contacts with nearly 50% of the C and F strands of CD147. On the other domain of CD147, parts of the A and G strands forms a complementary surface to the cyclic highly solvent-accessible segment C480-C488 of the spike RBD.

In the predicted complex, a row of hydrophilic residues along the interface formed strong and intricate polar contacts ( Figure 3 ). The most important interacting residues from both sides, together with the distances between each pair, are listed in Table 2 . Strong polar contacts were found to involve the spike RBD residues Gln493 and Ser494 with the side chain of Thr135 in CD147. Also, the side chain nitrogen of the spike Gln493 makes a 3.00 Å H bond with the hydroxy group of Ser190 in CD147. In addition, the backbone carbonyl of spike Glu484 formed a strong H bond (2.5 Å) with the backbone NH of Gln100 in CD147. Moreover, the side chain NH 2 of Asn481 on the cyclic segment of the spike b1'/b2' loop formed a bidentate H-bonds with the backbone of Phe27 and side chain OH of Thr28 in CD147. Another set of adjacent hydrophilic residues located on the spike RBD external domain (Gln498, Thr500, Asn501, Gly502 and Tyr505) were found to contribute a network of H-bonds with CD147 amino acids Trp137, Thr140, Ser163 and Thr188, which further augments the two proteins interaction. Furthermore, two strong salt bridges were observed to stabilize the binding interface at both ends of the spike F486- Interface area, calculated as difference in total accessible surface areas of isolated and interfacing structures divided by two. c D i G P-value indicates the P-value of the observed solvation free energy gain. The P-value measures the probability of getting a lower than observed D i G, when the interface atoms are picked randomly from the protein surface. p < 0.5 indicates interfaces with surprising (higher than would-be-average for given structures) hydrophobicity, implying that the interface surface can be interaction-specific.

G502 segment. The first one was noticed between Arg403 and Asp136 of the CD147 and the other was a stronger interaction (2.7 Å) emerging from spike Glu484 towards CD147 Lys191. The latter was acting as a centrally located anchor point that might have a key role in stabilizing the complex. It is worth noting that the same residue, Lys191, is an important player in the interaction of CD147 with the P. falciparum protein (PfRH5) (Wright et al., 2014) . Notably, fewer hydrophobic interactions were noticed at the interface as shown in Table 2 . The most important of these interactions was the one between the spike Tyr489 and His102 on the linker between the two domains of CD147.

In addition, to further validate the proposed binding mode, we have docked the spike glycoprotein of SARS-CoV, the causative agent of the 2003 epidemic (PDB ID: 5XLR) to the target protein CD147. This protein is known to bind with the human cell CD147 to facilitate entry of the virus responsible for the 2003 outbreak. We used HADDOCK server and the same parameters exploited for docking of the SARS-CoV-2 spike. To our delight, this protein adopted a very similar binding pose to that of the SARS-CoV-2 predicted above, with a slight shift towards the C-domain of the CD147 ( Figure S2 , Supp. Info.). Also, CD147 protein made polar contacts to the SARS-CoV protein using the same residues, namely, Gly100, His102, and Lys191. This gives added validity to our proposed binding pose of the SARS-CoV-2 spike and CD147.

The complex was then subjected to a Molecular Dynamics (MD) simulation study in order to check the stability of the proposed binding mode and to study the time-dependent behavior of this system. In absence of 3D structure information, MD simulations is especially useful for studying proteinprotein complexes as it can help the molecules "explore" their conformation space more efficiently than using the static "image" provided by the molecular mechanics energy minimization (Karplus & Petsko, 1990; Zhao et al., 2013) . The system was solvated in a cubic SPC water box and simulated for 100 ns in GROMACS software package under NPT ensemble using Nos e-Hoover thermostat (Van Der Spoel et al., 2005) . The final frame was extracted, minimized to a gradient of 0.001 kcal.Å À2 .mol À1 in MOE software, and was further analyzed for key structural changes (Vilar et al., 2008) . To our delight, the overall fold of both interacting proteins was preserved as well as most of the binding interface with an RMSD of just 1.25 Å. The spike RBD has moved slightly towards the CD147 N-terminal domain, while most of the polar interactions were retained or replaced by similar interactions with nearby residues (Figure 4(A, B) ). The backbone RMSD fluctuation of the whole complex as well as the interface residues were monitored during the MD simulation to check the stability of the system (Figure 4(C) ). Both proteins were stable with some RMSD fluctuations that stabilized after 40 ns. To our delight, the residues of CD147 at the interface attained a significant stability over time and reached 1 Å after 60 ns. Distances of the most important polar contacts at the interface were calculated throughout the simulation and represented as a heatmap in Figure 4(D) . It was noticed that the salt bridge between Arg403 of the spike and Asp136 of CD147 was stable during the simulation with only minor fluctuation for a short time between 20 and 36 ns. This ionic bond could be critical for anchoring the b-strand B of the CD147 C-terminal domain to the binding grove created by the spike RBD. Similarly, the H bond between Gln498 and Thr188 at the end of stand F of CD147, although weak, could contribute to some extent to the interaction of this strand to the surface of the spike. On the other side, at the beginning of the spike F486-G502 segment, the H bond of Asn481 to Thr28 was replaced by a stronger and more stable one with Thr44. The most prominent difference was the loss of the characteristic salt bridge between Glu484 and Lys191. In the starting complex, this interaction was centrally located below the loop connecting the two domains of CD147 and we initially thought it would be essential for stabilizing the interaction. As a compensation, both residues Glu484 and Lys191, upon liberation from the ionic bond, were able to form H bonds with Gln100 and Tyr449, respectively. Finally, we noticed that Tyr140 adopted a swing movement to bind with the backbone NH of Gly502 instead of Thr500, thus, stabilizing the strand D of CD147 to the edge of the spike F486-G502 segment.

The residues proposed to be critical for the Spike-CD147 interaction were mutated to alanine to validate their important role in the interaction. These residues are: Arg403 and Asn481 in spike, and Tyr140 in CD147. The mutated proteins were docked using HADDOCK server with the same parameters used for docking of the wild-type proteins. As expected, the number of residues at the interface and the interface surface area were lower than those of the wild-type proteins (Table S4 , Supp. Data). Also, the MM/GBSA score, as predicted by the Hawkdock server, was À60.64 which is lower than the scores of the wild-type protein complexes generated by any of the docking servers we used. Overall, the mutant proteins bind in a similar fashion to the wild-type one with an RMSD of 1.332 Å. However, the critical salt bridges and H bonds formed by the mutated residues were replaced by weaker contacts involving Phe486, Tyr489, Gln493, and Asn501 from spike with Val160, Asp136, Ser109, and Gln100 from CD147, respectively ( Figure S3 , Supp. Data).

In addition, the complex with the mutated residues was energy minimized and then subjected to MD simulation for 30 ns following the same protocol used for the wild-type complex. The purpose of this simulation is to compare the time-dependent behavior of both proteins and validate our proposed binding mode. It was noticed that the H bonds Tyr489-Asp136 and Gln493-Ser109 were maintained throughout the simulation. However, the complex was significantly less stable compared to the wild-type complex. This was evident from the higher fluctuation of the backbone RMSD of the whole complex and, especially, the interface residues ( Figure S4 , Supp. Data). These findings highlight the role of the mutated polar residues, Arg403 and Asn481 in spike, and Tyr140 in CD147, for the interaction of these proteins and give more confidence to our proposed binding mode.

To get more insight into the nature of interaction between the two proteins and detailed information about individual ligand contribution, we have performed binding free energy calculation (Cavasotto, 2020) . We have used the MD-based method Molecular Mechanics Poisson Boltzmann Surface Area (MM/PBSA) as implemented in the g_mmpbsa tool. This technique can achieve similar accuracy to the Free Energy Perturbation (FEP) methods but at a smaller computational cost (Kumari et al., 2014) . The calculation was performed with the single trajectory approach using 100 frames saved at regular intervals of the complex MD simulation. According to the MM/PBSA method, the binding free energy of the interacting proteins is calculated using the following equations:

where G complex is the total free energy of the complex and G protein is the total free energy of the isolated proteins. hE MM i is the vacuum molecular mechanics potential energy. TS refers to the entropic contribution to the free energy where T and S denote the temperature and entropy, respectively. The last term hG solvation i is the free energy of solvation and is composed of the polar part G polar and the nonpolar part G non-polar ; the former could be estimated by solving the Poisson À Boltzmann (PB) equation and the latter could be calculated from the solvent-accessible surface area (SASA). The term c is a coefficient related to surface tension of the solvent and b is fitting constant (Kumari et al., 2014) . The results of the binding free energy calculation are shown in Table 3 . In accordance with our initial prediction, the binding energy is dominated by the electrostatic term with a contribution of 66% compared to 30% of the van der Waals interaction. Therefore, the attractive electrostatic term greatly exceeded the repulsive term and was the main driving force for the binding. The MM/PBSA method also returned a Perresidue binding-free energy decomposition ( Figure S5 , Supp. Data). This helps in identifying the key residues involved in the interaction. In our case, a contribution of more than 5 kJ/ mol to the binding-free energy of the complex was considered significant. As expected, most of the residues participated in the binding of the initial and MD-refined complexes showed a significant contribution. Particularly, the calculation highlighted the role of the spike Arg403, Tyr449, Asn481, Val483, and Tyr505 in the attachment to the CD147 receptor. The overall proposed binding mode of the host cell CD147 and the SARS-CoV-2 spike is illustrated in Figure 5 .

Lymphopenia was found to be associated with the severity of SARS-CoV-2 infection in many clinical reports (Guan et occurrence of lymphopenia corresponds to a nearly 3-fold increased risk of severe disease , supporting its potential usage as a clinical sign for serious morbidity . Although the cause of lymphopenia during viral infection still needs further investigations, several studies have postulated potential mechanisms for the evident lymphopenia. These include direct viral invasion into T lymphocytes resulting in cell death , indirect upregulation of P53 apoptotic pathways in T lymphocytes by excessive secretion of pro-inflammatory chemokines and cytokines (IL10, IL7, IL6, IL2, TNFa, CXCL10/IP-10, CCL3/MIP-1A CCL2/MCP-1, and CCL4/MIP1B) Liao et al., 2002; Xiong et al., 2020) , and hyperlactic-acidemia-induced apoptosis (D ıaz et al., 2018) . Notably, T lymphocytes barely express ACE-2 receptor (Hamming et al., 2004 ) -the currently accepted point of entry of SARS-CoV-2 . CD147 may thus present an alternative point for viral spike protein entry. Indeed, CD147 was reported to mediate the entry of human immunodeficiency virus (HIV) into T lymphocytes (Pushkarsky et al., 2001) . It was also reported to mediate entry of cytomegalovirus (CMV) into endothelial and epithelial cells (Vanarsdall et al., 2018) , and Plasmodium falciparum into red blood cells (Zhang et al., 2018) . CD147, also termed as EMMPRIN or basigin, is a pleiotropic molecule with various effects (Grass & Toole, 2016; Iacono et al., 2007) . Given the emerging use of bioinformatics in biological research (Karim et al., 2020; Mehmood et al., 2014) , it is worth mentioning that, through pathway analyses, CD147 was found to be a significant cell-surface protein contributing to chemoresistance and cell survival in cancerstem-cell-like cells (Kang et al., 2013) , and the invasion of various cancers (Li et al., 2015; Liang et al., 2016; Zheng & Gong, 2017) . Additionally, it was found to be highly expressed in activated regulatory T lymphocytes, enabling it to be a potential marker for distinguishing cytokine-producing activated T lymphocytes from resting T lymphocytes (Landskron & Task en, 2013; Solstad et al., 2011) . The structure of transmembrane CD147 receptor consists of three domains: extracellular, transmembrane, and a short cytoplasmic tail (Muramatsu, 2016) . Those domains facilitate its interaction with several proteins, such as monocarboxylate transporters (MCTs) and cyclophilins (Cyps), especially cyclophilin A, to modulate specific functions (Kirk et al., 2000; Yurchenko et al., 2010) . CypA is a ubiquitously expressed chaperon protein that has a high affinity to bind to CD147 receptor for enhancing viral invasion and inflammation (Chen et al., 2005; Yurchenko et al., 2010) ; a study in 2013 reported that blocking of CypA inhibits SARS-CoV infection and replication (Tanaka et al., 2013) . CD147 acts as the major chaperone protein of MCT 1, MCT 3 and MCT 4 proton transporters, that control the export of lactic acid molecules outside the cell (Kirk et al., 2000; Philp et al., 2003) . Activated T lymphocytes have been reported to undergo increased glycolysis with increased lactate production (Cham & Gajewski, 2005) , as they need high energy supplies for proliferation and producing cytokines (Frauwirth & Thompson, 2004) . The expression of MCTs for lactate export helps T lymphocytes to cope with the increased lactic acid levels (Merezhinskaya et al., 2004) , since the accumulation of lactic acid induces T cells apoptosis (D ıaz et al., 2018) . These data provide an insight into the possible impact of SARS-CoV-2 infection on the balance of CD147-MCT interaction, and direct a need for further research on this role.

Angiotensin-converting enzyme 2, ACE-2 receptor, has been proved to be the entry point of SARS-CoV-2 into the lungs and other types of human cells. However, T lymphocytes are rarely expressing this receptor. Several recent studies have suggested that CD147 is a probable route for SARS-CoV-2 invasion into certain host cells. These findings encouraged us to design the current study in which we employed docking and MD simulation to investigate this possibility. The results were also compared with the molecular modeling data of the mutant proteins interaction as well as the interaction of CD147 with the spike of the SARS-CoV, the causative agent of the 2003 outbreak.

Our model shows interaction between CD147 C-terminal domain with SARS-CoV-2 spike external subdomain in the grove between the short antiparallel b strands, b1' and b2', and the small helix a1'. We confirmed this proposition using MD simulation and binding free energy calculation. This model provides an appreciation for targeting CD147 for adjuvant drug therapy strategies, to improve the prognosis of COVID-19 pateint. Our data is supported by the results of using Meplazumab, an anti-CD147 humanized antibody, for viral clearance and lymphocyte count normalization (Bian, 2020) . In addition, of proposing azithromycin as a therapeutic strategy to treat P. falciparum via interfacing in PfRH5/ CD147 interaction (Muralidharan & Striepen, 2015; Wilson et al., 2015) . The mechanisms proposed in our study encourage research on CD147 involvement in SARS-CoV-2 in direct invasion to T lymphocytes and associated lymphopenia.

No potential conflict of interest was reported by the author(s). 

Meplazumab treats COVID-19 pneumonia: An openlabelled, concurrent controlled add-on clinical trial

The coronavirus spike protein is a class I virus fusion protein: Structural and functional characterization of the fusion core complex

The species severe acute respiratory syndrome-related coronavirus: Classifying 2019-nCoV and naming it SARS-CoV-2

Binding free energy calculation using quantum mechanics aimed for drug lead optimization

Glucose availability regulates IFNgamma production and p70S6 kinase activation in CD8þ effector T cells

Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: A descriptive study

Function of HAb18G/ CD147 in invasion of host cells by severe acute respiratory syndrome coronavirus

Middle East respiratory syndrome coronavirus efficiently infects human primary T lymphocytes and activates the extrinsic and intrinsic apoptosis pathways

Basigin is a receptor essential for erythrocyte invasion by Plasmodium falciparum

Colocalization of galectin-3 with CD147 is associated with increased gelatinolytic activity in ulcerating human corneas

CPORT: A consensus interface predictor and its performance in prediction-driven docking with HADDOCK

Unravelling the interplay between extracellular acidosis and immune cells

HADDOCK: A proteinprotein docking approach based on biochemical or biophysical information

Inhibitory effect of tumor cell-derived lactic acid on human T cells

Regulation of T lymphocyte metabolism

Dynamic attractor for the Berendsen thermostat an the slow dynamics of biomacromolecules

Severe acute respiratory syndrome-related coronavirus-The species and its viruses, a statement of the Coronavirus Study Group

How, with whom and when: An overview of CD147-mediated regulatory networks influencing matrix metalloproteinase activity

Clinical characteristics of coronavirus disease 2019 in China

Cryo-electron microscopy structures of the SARS-CoV spike glycoprotein reveal a prerequisite conformational state for receptor binding

Tissue distribution of ACE2 protein, the functional receptor for SARS coronavirus. A first step in understanding SARS pathogenesis

Combined receptor-based and ligandbased approach to delineate the mode of binding of guaianolideendoperoxides to PfATP6

New insights into the binding mode of melanin concentrating hormone receptor-1 antagonists: Homology modeling and explicit membrane molecular dynamics simulation study

SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor

Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations

VMD: Visual molecular dynamics

CD147 immunoglobulin superfamily receptor function and role in pathology

COVID-19 spike-host cell receptor GRP78 binding site prediction

CD147 as a key mediator of the spleen inflammatory response in mice after focal cerebral ischemia

Proteomic analysis reveals that CD147/ EMMPRIN confers chemoresistance in cancer stem cell-like cells

A multi-omics analysis of Bone Morphogenetic Protein 5 (BMP5) mRNA expression and clinical prognostic outcomes in different cancers using bioinformatics approaches

Molecular dynamics simulations in biology

Human leukocyte activation antigen M6, a member of the Ig superfamily, is the species homologue of rat OX-47, mouse basigin, and chicken HT7 molecule

The E-selectin ligand basigin/ CD147 is responsible for neutrophil recruitment in renal ischemia/ reperfusion

CD147 is tightly associated with lactate transporters MCT1 and MCT4 and facilitates their cell surface expression

T cell activation-associated epitopes of CD147 in regulation of the T cell response, and their definition by antibody affinity and antigen density

Inference of macromolecular assemblies from crystalline state

g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

Severe lymphopenia in hospitalized patients with influenza virus infection as a marker of a poor outcome

Impact of severe hematological abnormalities in the outcome of hospitalized patients with influenza virus infection

CD147 in regulatory T cells

Structure, function, and evolution of coronavirus spike proteins

CD147 reprograms fatty acid metabolism in hepatocellular carcinoma cells through Akt/mTOR/SREBP1c and P38/PPARa pathways

Expression of the SARS-CoV-2 cell receptor gene ACE2 in a wide variety of human tissues

E2F1 promotes tumor cell invasion and migration through regulating CD147 in prostate cancer

IL-19 induces production of IL-6 and TNF-alpha and results in cell apoptosis through TNF-alpha

Analysis of factors associated with disease outcomes in hospitalized patients with 2019 novel coronavirus disease

Clinical and biochemical indexes from 2019-nCoV infected patients linked to viral loads and lung injury

Outbreak of pneumonia of unknown etiology in Wuhan China: The mystery and the miracle

Genomic characterisation and epidemiology of 2019 novel coronavirus: Implications for virus origins and receptor binding

Use of bioinformatics tools in different spheres of life sciences

Presence and localization of three lactic acid transporters (MCT1,À 2, andÀ 4) in separated human granulocytes, lymphocytes, and monocytes

Teaching old drugs new tricks to stop malaria invasion in its tracks

Basigin (CD147), a multifunctional transmembrane glycoprotein with various binding partners

Loss of MCT1, MCT3, and MCT4 expression in the retinal pigment epithelium and neural retina of the 5A11/basigin-null mouse

ZDOCK server: Interactive docking prediction of protein-protein complexes and symmetric multimers

CD147 facilitates HIV-1 infection by interacting with virus-associated cyclophilin A

World Health Organization declares global emergency: A review of the

CD147 (Basigin/Emmprin) identifies FoxP3 þ CD45RO þ CTLA4þ-activated human regulatory T cells

Cryo-EM structure of the SARS coronavirus spike glycoprotein in complex with its host cell receptor ACE2

SARS-CoV-2 entry factors are highly expressed in nasal epithelial cells together with innate immune genes

Lymphopenia predicts disease severity of COVID-19: A descriptive and predictive study

Suppression of coronavirus replication by cyclophilin inhibitors

CASTp 3.0: Computed atlas of surface topography of proteins

An approach to developing a force field for molecular simulation of martensitic phase transitions between phases with subtle differences in energy and structure

Extracellular matrix metalloproteinase inducer EMMPRIN (CD147) in cardiovascular disease

GROMACS: Fast, flexible, and free

CD147 promotes entry of pentamer-expressing human cytomegalovirus into epithelial and endothelial cells

Medicinal chemistry and the molecular operating environment (MOE): Application of QSAR and molecular docking to drug discovery

SARS-CoV-2 invades host cells via a novel route: CD147-spike protein

Structural and functional basis of SARS-CoV-2 entry by using human ACE2

Clinical features of 69 cases with coronavirus disease 2019 in Wuhan, China

HawkDock: A web server to predict and analyze the protein-protein complex based on computational docking and MM/GBSA

Macrolides rapidly inhibit red blood cell invasion by the human malaria parasite, Plasmodium falciparum

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Structure of malaria invasion protein RH5 with erythrocyte basigin and blocking antibodies

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients

High expression of ACE2 receptor of 2019-nCoV on the epithelial cells of oral mucosa

Cyclophilin-CD147 interactions: A new target for anti-inflammatory therapeutics

Mortality of COVID-19 is associated with cellular immune function compared to immune function in Chinese Han population

Clinical characteristics of 140 patients infected with SARS-CoV-2 in Wuhan

Disrupting CD147-RAP2 interaction abrogates erythrocyte invasion by Plasmodium falciparum

Mature HIV-1 capsid structure by cryo-electron microscopy and all-atom molecular dynamics

Lymphopenia is associated with severe coronavirus disease 2019 (COVID-19) infections: A systemic review and meta-analysis

CD147 expression was positively linked to aggressiveness and worse prognosis of gastric cancer: A meta and bioinformatics analysis

Functional exhaustion of antiviral lymphocytes in COVID-19 patients

Improving the performance of the MM/PBSA and MM/GBSA methods in recognizing the native structure of the Bcl-2 family using the interaction entropy method

A pneumonia outbreak associated with a new coronavirus of probable bat origin

The authors would like to thank the Bibliotheca Alexandrina, Alexandria, Egypt for providing access to the High-Performance Computing Cluster which was used for the molecular dynamics simulation part of the study, ASRT grant #7301 for partial funding, and internal funding ZC003-2019 from Zewail City for Science and Technology.