Predicting a Kind of Unusual Multiple-States Dimerization-Modes Transformation in Protein PD-L1 System by Computational Investigation and a Generalized Rate Theory

The new cancer immunotherapy has been carried out with an almost messianic zeal, but its molecular basis remains unclear due to the complexity of programmed death ligand 1 (PD-L1) dimerization. In this study, a new and integral multiple dimerization-modes transformation process of PD-L1s (with a new PD-L1 dimerization mode and a new transformation path discovered) and the corresponding mechanism are predicted using theoretical and computational methods. The results of the state analysis show that 5 stable binding states exist in system. A generalized inter-state transformation rate (GITR) theory is also proposed in such multiple-states self-assembly system to explore the kinetic characteristics of inter-state transformation. A “drug insertion” path was identified as the dominant path of the PD-L1 dimerization-modes transformation. Above results can provide supports for both the relative drug design and other multiple-states self-assembly system from the theoretical chemistry perspective.


INTRODUCTION
The search for a cure for cancer has been one of the great pursuits of modern science. Cancer immunotherapy has been called the third revolution in cancer therapy (Miller and Sadelain, 2015), bringing hope for a cure for many types of cancer (Mahoney et al., 2015;Hoos, 2016;Hu and Li, 2020). One of the main mechanisms of cancer immunotherapy is to block the binding between the receptor, known as programmed death protein 1 (PD-1), on the surface of immune cells and the ligand, known as programmed death ligand 1 (PD-L1), on the surface of cancer cells to prevent cancer cells immune evasion, so as to promote the recognition and the targeted killing of cancer cells by the body's immune cells (Ishida et al., 1992;Sharma and Allison, 2015).
The PD-1/PD-L1 target of cancer immunotherapy has already exhibited excellent clinical pharmacological effects. For example, it has been a great success to block the binding of PD-1/ PD-L1 with monoclonal antibody in clinical application Topalian et al., 2012). However, the antibody drugs have some inherent disadvantages such as high price, short half-life, poor penetrability and immunogenicity (Chames et al., 2009). Fortunately, these disadvantages can be overcome by the development of small-molecule PD-1/PD-L1 blocking drugs. It is a pity that the development of small-molecule drugs that block the PD-1/PD-L1 pathway is relatively lagging currently, which is mainly due to the lack of in-depth understanding of the interaction mechanism between the small-molecule drugs and the PD-L1s.
The Bristol-Myers-Squibb (BMS) team has made a breakthrough in developing a class of effective small-molecule PD-1/PD-L1 blocking drugs based on the (2-methyl-3biphenylyl) methanol scaffold (Chupak et al., 2017;Chupak and Zheng, 2018), and the mechanism of interaction between the drugs and the proteins had not been reported. Holak et al. solved several holo form (involving small-molecule blocking drugs) sandwich PD-L1 dimeric crystal structures (PDB ID: 5J89, 5J8O, 5N2D, 5NIU, 5N2F, 6R3K, 6RPG) using X-ray (Zak et al., 2016;Guzik et al., 2017;Skalniak et al., 2017;Basu et al., 2019). These holo form dimerization structures in the solid phase crystals provide a prerequisite for relevant studies and greatly contribute to the understanding of the action mechanism of such small-molecule PD-1/PD-L1 blockers. More interestingly, to apo form PD-L1, there are several X-ray diffraction crystal structures (PDB ID: 4Z18, 5JDR, 3FN3) and EGS cross-linking molecular weight evidence indicating that PD-L1s may also be apo form dimer in the absence of drugs (Chen et al., 2010;Zhang et al., 2017). Even so, PD-L1s have different binding interfaces and nearly opposite relative orientation between holo form dimerization mode and apo form dimerization mode.
The above results indicate that PD-L1s may have both holo form dimers and apo form dimers in solution. Together, these results provide complete prerequisite information to explore the integral mechanism of drug involving PD-L1 dimerization, for that it is difficult to recruit macromolecule PD-L1s aggregating entirely by small molecules themselves in physical images. However, currently, the potential transformations between the apo form dimerization modes and the holo form dimerization modes remain ambiguous, as none of the existing studies has adequately considered the forming of holo dimers in relation to the existence of the apo dimers.
To clarify the remaining research space, we address the following questions: To begin with, in the aspect of stable states of the system, there are three questions below. Whether the two kinds of dimerization modes (holo form and apo form) solved using X-ray in solid phase can be stable in solution? Are there any other stable binding states in the solution? What are the stability factors for dimerization modes that can exist in solution? Further, in the aspect of inter-state transitions among the stable states, there are also three more questions. What are the possible self-assembly paths that can happen in kinetics? How to determine the dominant path in this multiple-paths selfassembly system? Ultimately, does this class of small-molecule drugs act as inducers or stabilizers in the formation of the holo form dimerization mode? Combined computations of multi-scale simulations with theoretical deduction, the above questions could be solved, which will, in turn, bridge the gaps within the current research.
Moreover, for the kinetic characteristics in the abstract aspect, such as the dominant path of inter-state transformation in this complex self-assembly network, we proposed a generalized interstate transformation rate (GITR) theory in a different perspective from the numerical calculation. This work provides additional insights into other similar complex multiple-states self-assembly system.

Starting Structures in Molecular Dynamics Simulations
The initial structure of apo form dimerization structure 1 was obtained from Protein Data Bank (PDB ID: 4Z18). The two chains of 4Z18 were cut off at the peptide bonds between Tyr134 and Asn135, with the relative position and orientation of the two chains kept unchanged. And then, the upper IgV domain (D1) and the lower IgC domain (D2) of PD-L1 from the initial structure 4Z18 were obtained and renamed as 4Z18 U and 4Z18 L respectively (Chen et al., 2010). Correspondingly, the initial structure of holo form dimerization structure was obtained from Protein Data Bank (PDB ID: 5J89) (Zak et al., 2016). Similarly, the starting structures of IS L , IS R and IS PP were obtained by removing A chain, B chain and drug molecule BMS-202 from the structure of 5J89 respectively, keeping the relative position and orientation of the rest parts unchanged. The H++ online server (Gordon et al., 2005) was used to calculate the protonation states of the ionizable residues at pH 7.0. The partial charges of drug molecule BMS-202 were calculated by AM1-BCC method (Jakalian et al., 2000). The field parameters of BMS-202 were obtained by using Antechamber Suite with the General AMBER Force Field (GAFF) (Wang et al., 2004). Similarly, the field parameters of PD-L1 were obtained by using t-LEaP module of AMBER 16 package with the ff14SB force field (Maier et al., 2015;Case et al., 2016). Moreover, the missing atoms in the proteins were added by t-LEaP module. In order to make the system electrically neutral, Na + or Cl − models were added as counter ions, with their three-dimensional positions were computed using Coulombic potential grid algorithm (Case et al., 2016). Finally, each system was solvated using the TIP3P water model (Jorgensen et al., 1983) in a truncated octahedron box with a 12.0 Å distance around the solute.

Molecular Dynamics Simulations
The MD simulations (Alder and Wainwright, 1959;Mccammon et al., 1976;Tuckerman and Martyna, 2000) were carried out using the AMBER 16 software package (Case et al., 2016). Firstly, 20,000 steps of energy minimization (10,000 steps of steepest decent followed by 10,000 steps of conjugate gradient) were carried out with protein and drug constrained (500 kcal mol −1 Å −2 ). Subsequently, 20,000 steps of energy minimization were repeated without any constrain. Next, each system was heated linearly from 0 to 310 K over a period of 310 ps with 10.0 kcal mol −1 Å −2 restrain on the solute and then another 200 ps equilibrium simulation in the NVT ensemble was followed at 310 K with 5.0 kcal mol −1 Å −2 restrain on the solute. Finally, 150 ns MD simulation was performed for each system in the NPT ensemble to produce trajectory. A constant isotropic pressure of the system was maintained at 1 atm using the Berendsen barostat (Berendsen et al., 1984). Meanwhile, the temperature of each system was maintained at 310 K by coupling to a Langevin heatbath (Uberuaga et al., 2004) using a collision frequency of 1 ps −1 . Short range nonbonded interactions were cut off at 10.0 Å, while the long-range electrostatic interactions were handled using the particle mesh Ewald (PME) method (Darden et al., 1993). The time step of each MD simulation was set to 2 fs using SHAKE algorithm (Ryckaert et al., 1977) to restrict all covalent bonds involving hydrogen atoms.

Adaptive Steered Molecular Dynamics Simulations
The starting state of ASMD simulations (Ozer et al., 2010;Ozer et al., 2012) was the last frame in the 150 ns MD simulation of the crystal structure (PDB ID:5J89). The adaptive steering displacement vector was defined with the ALA-B121 Beta C atom in the PD-L1s dimer as the starting point and the C12 atom at the end of drug molecule BMS-202 as the terminal point. The length of the adaptive steering displacement vector was measured by VMD (Humphrey et al., 1996). Considering the length of the drug molecule, we pulled the drug along the displacement vector 20 Å with the magnitude of the vector increasing from 17.28 to 37.28 Å. The pulling speed v and the spring constant k were chosen as 2 Å ns −1 and 20 kcal mol −1 Å −2 respectively. The potential of mean force (PMF), representing the free energy change ΔG along the displacement vector or pulling coordinate, would be computed from ASMD trajectories using Jarzynski's equality (Ozer et al., 2010;Ozer et al., 2012): where the W s → t is the pulling work from the starting point to the terminal point in trajectories. In order to adapt our ASMD algorithm to the geometry of the sinuous channel between two PD-L1 monomers, we divided 20-Å-pulling coordinate into 10 stages to increase the calculation accuracy. At each stage, steering molecular dynamics simulations were run 20 times independently to enhance the sampling, with 500,000 steps for each trajectory and 0.002 ps for each step. Totally, for each system, 200 ns ASMD simulations were needed to pull the drug out of the channel. The ASMD simulations were executed by the AMBER package (Case et al., 2016).

Root Mean Square Deviation and Root Mean Square Fluctuations
Root mean square deviation measures the change of the conformation with the change of time, while root mean square fluctuations measure the fluctuations of atoms or residues during the simulations. They are defined below: Here, the r i (t) is the position vectors of atom i at time t and the r i ′ is the position vectors of atom in reference structure. Specifically, the N, Frames and R represent summing over each atom in the selected parts of the structure, each frame in the trajectories and each atom in the residue respectively. They were computed using the CPPTRAJ module of the AMBER package (Roe and Cheatham, 2013).

Dynamic Cross-Correlation Matrix
Firstly, all the conformations in the trajectories would be aligned against a reference structure to remove overall global translational and rotational motions by means of a least-square-fitting procedure using all backbone Cα atoms. Covariance matrix elements are denoted as lower-case c ij while dynamical crosscorrelation matrix elements are denoted as upper-case C ij . They are computed by Eqs 5, 6 as shown below (Hunenberger et al., 1995): Here, the angle brackets denote time averages. The r i (t) and r j (t) are the position vectors of atom i and atom j at time t respectively. Dynamical cross-correlation matrix element can be understood as the normalized covariance matrix element with the value from −1 to +1. It reflects the correlation of displacements between two atoms along a straight line. The DCCM were computed by the CPPTRAJ module of the AMBER package (Roe and Cheatham, 2013).

Binding Free Energy Calculation
Binding free energy G bind can be divided into G mmpbsa term (in which only solvation entropy is considered in the aspect of entropy effect, as defined below) and gas phase conformational entropy term TS, which could be calculated by MM-PBSA method and Normal Mode Analysis method (Miller et al., 2012) respectively: G mmpbsa can be divided into the molecular mechanics energy term E MM and the free energy change of solvation effect ΔG solv : Furthermore, the molecular mechanics energy term E MM can be divided into the internal energy of molecule E int , the van der Waals energy E vdw and the electrostatic energy E ele : Moreover, the internal energy of molecule E int can be divided into the bond energy term E bond , the angle energy term E angle , the torsion energy term E torsion : Correspondingly, the free energy change of solvation effect ΔG solv can be divided into the free energy change of polar solvation effect ΔG PB and the free energy change of nonpolar solvation effect ΔG SA : Note that the free energy change of polar solvation effect ΔG PB could be obtained by solving Poisson-Boltzmann equation (Fogolari et al., 1999): where ε( r) is the dielectric constant, Ø( r) is the electrostatic potential, ρ( r) is the solute charge, λ( r) is the Stern layer masking function, z i is the charge of ion type i, c i is the bulk number density of ion type i far from the solute, k B is the Boltzmann constant, and T is the temperature.
The free energy change of nonpolar solvation effect ΔG SA can be calculated from surface tension c and solvent accessible surface SASA Sitkoff et al., 1994. ΔG SA cSASA + β As a result, the change of binding free energy ΔG bind can be deduced from ΔG mmpbsa and the change of gas phase conformational entropy TΔS, as follows: where the superscripts A, B and AB denote receptor, ligand and complex of receptor and ligand respectively.

RESULTS AND DISCUSSION
The Stable States in the PD-L1s System These stable binding states can be divided into the initial states (apo PD-L1s dimerization), the final state (holo PD-L1s dimerization), and the intermediate states (stable binding modes in the transformation process of initial states to final state) according to the drug-regulated dimerization process.
Initial State 1: Apo PD-L1s Dimerization Mode in 4Z18 All the apo dimerization modes (binding modes) in the crystal structure (PDB ID: 4Z18, 5JDR, 3FN3) are same, and the crystal structure (PDB ID: 4Z18) has the best X-ray diffraction resolution (Chen et al., 2010;Zhang et al., 2017). Thus, the crystal structure (PDB ID: 4Z18) was used as the starting structure for the 150 ns all-atom molecular dynamics (MD) simulation (Alder and Wainwright, 1959;Mccammon et al., 1976;Tuckerman and Martyna, 2000) under the condition described as Molecular Dynamics Simulations methods. After visualizing the simulation trajectory using VMD program (Humphrey et al., 1996), it was found that the upper half (IgV domain) of the two PD-L1s had a slight oscillating motion of opening and closing, while the lower half (IgC domain) bound tightly all the time. That is to say, this dimerization mode did not dissociate in the simulation. The root mean square deviation (RMSD) shown in Supplementary Figure S1 was used to quantitatively characterize the stability of the double-chain system (black line: 4Z18_whole_dul) and the single-chain system (red line: 4Z18_whole_sin) respectively. The fluctuation of doublechain RMSD is larger than that of single-chain RMSD, which indicates that there is a certain inter-chain relative motion. This motion conforms to the protein dimer system itself characteristic. From the last 100 ns equilibrium trajectory, 2,500 frames were uniformly taken out and the molecular mechanics Poisson-Boltzmann surface area method (MM-PBSA) was used to calculate the binding free energy. The ΔG bind value of 4Z18 whole system is −23.63 ± 10.48kcal/mol (Supplementary Table S1). This indicates that the dimerization mode in crystal structure 4Z18 can be stable thermodynamically. Our results can systematically quantify and verify the inference of existence of apo form PD-L1 dimerization in the previous research (Chen et al., 2010).
The binding sites of PD-L1 with such drugs or the receptor PD-1 are all in the upper IgV domain, while the lower IgC Domain of PD-L1 is an extracellular supporting structure (Chen et al., 2010). Being inspired by the motion mode of "upper loose and lower tight" in our visualization, this 4Z18 dimerization mode was truncated to the upper part 4Z18 U and the lower part 4Z18 L with the relative position remained unchanged to further clarify the precise binding site of the 4Z18 dimerization mode. The 4Z18 U and the 4Z18 L were used as the initial structures for 150 ns MD simulation respectively. In 4Z18 U , the RMSD of double-chain and single-chain structure are shown in the black line and the red line in Supplementary Figure S2A respectively. The significant difference between the two lines reflects the obvious relative motion between the two chains in 4Z18 U . The same analysis on the stability of 4Z18 L shows that the 4Z18 L structure has high inter-chain stability, as shown in Supplementary Figure S2B. After visualization of the two trajectories, it was found that the two chains of 4Z18 U exhibited obvious opening-closing oscillation in the first 20 ns, and dissociated with significant relative displacement after 60 ns, while the two chains of 4Z18 L binding stably in the whole 150 ns trajectory.
Frontiers in Chemistry | www.frontiersin.org November 2021 | Volume 9 | Article 783444 In the thermodynamical aspect, the calculated values of ΔG bind of 4Z18 U and 4Z18 L are 22.71 ± 8.86kcal/mol (Supplementary Table S2) and −21.62 ± 9.66kcal/mol (Supplementary Table S3) respectively. As our systems are protein-protein interaction (PPI) (Pitre et al., 2008) systems, the standard derivations of the binding free energies ΔG bind may look a little bit larger than the common protein small molecular systems due to much larger interface conformational flexibility and fluctuation. This phenomenon can also be seen in other similar PPI systems (Nayeem et al., 2017;Shi et al., 2019).
The above results show that the binding site of the 4Z18 dimerization mode is mainly in the lower part of PD-L1, as the upper part has no binding effect in this relative orientation. These results have not been reported in relevant studies, which is maybe mainly because the current studies pay too much attention to the upper half structural domain of PD-L1, and use the truncation model which only preserves the upper part in relevant calculation studies (Zak et al., 2016;Guzik et al., 2017;Skalniak et al., 2017;Basu et al., 2019;Shi et al., 2019;Soremekun et al., 2019;Sun et al., 2019).

Final State: Holo PD-L1s Dimerization Mode in 5J89
The structure of BMS-202 is close to the common scaffold of this kind of (2-methyl-3-biphenylyl) methanol drugs and relatively simple and representative in them. Especially, experiments show that BMS-202 also has strong pharmacological activity (Chupak et al., 2017;Chupak and Zheng, 2018). As can be seen, an indepth study of the action mechanism of BMS-202 is particularly important to improve the efficacy of such kind of drugs in the future. For these reasons, we chose PDB 5J89 (which contains BMS-202) as the object in this study.
The crystal structure (PDB ID: 5J89) was used as the starting structure for the 150 ns all-atom MD simulation under the condition described in Molecular Dynamics Simulations methods. As shown in Supplementary Figures S3A,B, the 5J89 is a kind of sandwich structure, with the small-molecular drug BMS-202 entrapped between two PD-L1 monomers.

Global Analysis in Protein Level: Stability, Flexibility, and Energy
The double-chain RMSD of the 5J89 dimer is shown by the black line in Supplementary Figure S3A. Considering its relatively large fluctuations, the root mean square fluctuations (RMSF) of each residual in the 150 ns trajectory was calculated, and the results were shown in Supplementary Figure S3B. The RMSF shows that the fluctuations mainly come from the 10-residuelong disordered structure at the end of two chains, which comes from the imperfect of artificial expression of proteins and is not belong to homo sequence. Since it is far away from the binding sites of PD-L1 with drugs and the receptor PD-1, it has little influence on our study. Hence the RMSD of 5J89 homo part (Residue ID: 18-134) was calculated in the 150 ns trajectory, as shown by the red line in Supplementary Figure S3A. The results show that the 5J89 dimerization mode has no obvious inter-chain relative motion. The visualization also showed that this sandwich structure did not dissociate in the whole 150 ns trajectory. For this 5J89 model, when calculating the binding free energies, there are some considerations needed in terms of geometric symmetry. Because the two PD-L1 monomers are exactly the same sequence, under the ensemble statistics in the case of removal of the BMS-202 molecule, the spatial structures of the two PD-L1 conformations (on average) should be exactly the same, the combination structure of two PD-L1s in the conformational average should be C 2v symmetry. However, the average conformation of BMS-202 should have σ V symmetry.
When the X, Y and Z three-dimensional coordinate system is established for the conformational average structure of 5J89 system and let the XZ plane pass through the plane of the conformational average structure of BMS-202, namely the conformational average binding interface of the two PD-L1 monomers, then the symmetry of the conformational average structure of the two PD-L1 monomers satisfies C 2v C 2Z in this coordinate system, the symmetry of the conformational average structure of the BMS-202 satisfies σ V σ XZ , and their matrices in this coordinate system are as follows respectively: This means that only when the average conformational structure of BMS-202 also has σ yz symmetry, or when the conformational average of the combination structure of two PD-L1 monomers also has σ yz symmetry, the binding interface environment between the BMS-202 and the two PD-L1 monomers is completely the same. However, neither of these two conditions can be satisfied. Therefore, according to the above simple analysis of group theory, we need to study the two-binding interface between the BMS-202 and two PD-L1 monomers respectively.
For the reasons of symmetry mismatch (between C 2v and σ V ) mentioned above, we calculated ΔG bind of three kinds of stripping modes (5J89 A , 5J89 B and 5J89 C ) respectively, shown in Figure 1. The detail components of binding free energies are shown in Supplementary Tables S4, S5, S6 respectively.
The negative free energies of three kinds of stripping modes in final state 5J89 have shown the mutual affinity between the three self-assembly blocks (when they stochastic move to the neighbourhood among them), which indicates the feasibility of three-body concerted self-assembly. The results of stripping mode 5J89 A and stripping mode 5J89 B show that the symmetry mismatch indeed leads to the difference of binding free energy. The ΔG bind of 5J89 C indicates that BMS-202 has a stabilizing effect on PD-L1s dimerization. The ΔG bind of 5J89 A and 5J89 B are far higher than the ΔG bind of 5J89 C , which indicates that there is a strong protein-protein interaction between the two PD-L1 monomers in 5J89 mode. This is further confirmed in the later (Three Other Binding States section). Our binding free energy results provide details in stripping modes, and also consistent with other works about other similar but not the

Local Analysis in Residues Level: Energy Decomposition, Key Residue Groups, H-Bonds
To further quantify the local structure factors contributing to the stability of the 5J89 dimerization mode from the thermodynamic aspect, the MM-PBSA binding free energy decomposition was performed, and the results are shown in Figure 2. Residues that contribute significantly to the binding stability of the 5J89 dimerization model are labeled in Figure 2. Our energy decomposition results can quantify and verify the previous prediction of key residues in PD-L1 by the Holak. research group (Zak et al., 2016;Guzik et al., 2017). Combined with the simulation trajectory visualization, we can summarize them into the following four key residue groups, as shown in Figure 3.
Group 1: There are four tyrosine residues (Tyr-A56, Tyr-B56, Tyr-A123, Tyr-B123) on the binding surfaces of BMS-202 and PD-L1s, as shown in Figure 3A. We call them "tyrosine sucking disks." They can oscillate and contact the three hydrophobic aromatic rings of the BMS-202 in various conformations of the NPT ensemble sampling, which can provide stabilization effect by hydrophobic adsorption.
Group 2: Around the BMS-202, there are two isoleucine residues (Ile-A54, Ile-B54), two valine residues (Val-B68, Val-B76), and two methionine residues (Met-A115, Met-B115) branching out six long hydrophobic side chains, as shown in Figure 3B. We call them "hydrophobic shrimp feet," since a shrimp has three pairs of jaw feet in its mouth. They can oscillate and contact the hydrophobic skeleton of BMS-202 in various conformations of the NPT ensemble sampling, which can provide a hydrophobic chelating stabilization effect.
Group 3: There is a pair of alanine residues (Ala-A121, Ala-B121) clipping the two benzene rings of BMS-202 by σ-π interactions, as shown in Figure 3C. We call them "alanine clips." Although this structure is small, the steric effect of its hydrophobic geometric complementation plays an important role in preventing BMS-202 shedding.
Group 4: There are three pairs of salt-bridges which anchor to each other by static electricity attraction of opposite charges, that are (−)Glu-A58/(+)Arg-B125, (−)Asp-A61/(+)Arg-B113, (+) Arg-A125/(−)Glu-B58, as shown in Figure 3D. The importance of electrostatic effect in this system is not only reflected in the high contribution of binding free energy decomposition, but also in the fact that electrostatic force is a long-range force, which is extremely important in the selfassembly remote recognition stage.
To further reveal the electrostatic driving effect, APBS program (Baker et al., 2001) was used to calculate the electrostatic potential on the molecular surface of full-length PD-L1 protein, as shown in Figure 4. The density functional FIGURE 1 | Three kinds of stripping modes of ΔG bind calculation in 5J89 system which are denoted as 5J89 A (A), 5J89 B (B) and 5J89 C (C) respectively. The gray surfaces represent the stripping interfaces in binding free energy calculations. B3LYP was used to calculate the surface electrostatic potential of BMS-202 at the base group level of 6-311+G (2d,p) by Gaussian 09 package (Andersson and Uvdal, 2005;Frisch et al., 2009). The calculation results of the electrostatic potential of the molecular surfaces show that the C 2v symmetry generated by the relative orientation of PD-L1s and the special charge distribution on the surfaces of the PD-L1 upper IgV domains together lead to the electrostatic opposite-matching exactly. Moreover, the characteristic of "alternating electropositive and electronegative of acetylethylenediamine" in the tail of BMS-202 also matches the electrical property on the surface of PD-L1s. This suggests that enhancing the alternating electrostatic potential difference in the tail structure of drug molecules may be more effective in the remote recognition stage of electrostatic driving.
In order to investigate intermolecular H-bonds across binding interface, we conducted statistical analysis of intermolecular H-bonds in 15,000 frames of conformations in our 150 ns trajectory using geometric criteria (Jones et al., 2014) that distance is less than 3.5 Å and angle is greater than 120°and less than 180°. The statistical results are shown in Supplementary Table S7.
The total occupancy of intermolecular H-bond is 11.47, which means that any frame of conformation taken from the NPT ensemble will have 11.47 intermolecular H-bonds on average, which is very considerable in the contribution of intermolecular binding free energy. It is worth mentioning that the residues with high occupancy in the H-bond statistics also contribute highly to the decomposition of binding free energy.

Three Other Binding States
The path of a self-assembly system is not single in many cases (Hagan and Chandler, 2006;Grzybowski et al., 2017;Michaels et al., 2018), and therefore the feasibility of other self-assembly path is still worthy to investigate. If the self-assembly process of final state 5J89 is a step by step process, then all the three possible formation processes can be as follows: When we choose a view pose that the acetamide tail of BMS-202 pointing towards us, 1) the PD-L1 binds with BMS-202 as left monomer firstly, 2) the PD-L1 binds with BMS-202 as right monomer firstly, 3) two PD-L1 monomers bind firstly. The three binding states are named as IS L , IS R and IS PP respectively and shown in Supplementary Figure S4. Similar to Global Analysis in Protein Level: Stability, Flexibility, and Energy section, symmetry mismatch (the C 2v group and σ v group) results in different relative binding orientation of the BMS-202 in IS L and IS R , which also needs to be studied separately. IS L , IS R and IS PP structures were simulated for 150 ns all-atom MD and binding free energies were calculated under the same conditions as before. The results are shown in Supplementary Figures S4A,B,C respectively. The detail components of binding free energies are shown in Supplementary Tables S8, S9, S10 respectively. In the visualization of the MD trajectories, none of the three binding states dissociated. The above results show that all the three binding states can exist stably, hence the step-by-step self-assembly is feasible. The three stable binding states are discussed in detail as follows.
IS R Is More Stable Than IS L as Asymmetry IS R has −0.14kcal/mol binding free energy superior to IS L , which is also consistent with the larger binding free energy when the right PD-L1 monomer stripping from 5J89 as calculated in Final State: Holo PD-L1s Dimerization Mode in 5J89 section. Combined with the 150 ns trajectories of IS L and IS R and the results of binding free energy decomposition, it can be found that the hydrophobic side chains of Val-B68 (−1.55kcal/mol) and Val-B76 (−1.19kcal/mol) in IS R are able to chelate the hydrophobic tail of BMS-202. As is shown in Supplementary Figure S5. It is easy to know from the rotational symmetry operation of C 2v that the corresponding Val-A68 and Val-A76 in IS L locate on the distal side of the benzene ring head of BMS-202, and they cannot provide the same stabilization as IS R in geometrically and hydrophobic environments. That is to say, asymmetry is what makes the difference.

A New Apo Form Dimerization Mode IS PP
We discovered a new kind of stable apo form dimerization mode (IS PP ), which is different from the apo form dimerization 4Z18 reported by previous research (Chen et al., 2010). The result of binding free energy shows that IS PP can exist as a stable dimerization mode in thermodynamic aspect, and compared FIGURE 4 | The surface electrostatic potential of PD-L1 and drug molecule BMS-202. The relative orientation of PD-L1s in 5J89 system satisfies the C 2v symmetry, and it will satisfy the translation symmetry if we rotate the left monomer 90°clockwise and the right monomer 90°counterclockwise. For ease of viewing, the electrostatic opposite-matching regions of the PD-L1s are linked in pairs by black lines. The surface electrostatic potential and the 2D structure with the atomic numbers of the drug BMS-202 are shown respectively.  (Miller et al., 2012). As shown in Table 1, the ΔE nonpolar in 4Z18 dimer is high and plays a stabilization effect, and the ΔE polar in 4Z18 plays a repulsion effect. While in IS PP the ΔE polar plays a stabilization effect. Therefore, experimental crystallization conditions with high ionic strength may enhance the dimerization stability of 4Z18, but weaken the stability of IS PP , making it difficult to find IS PP in the crystallization state.

BMS-202 Can Stabilize the IS PP
Comparing the RMSD analysis of IS PP with 5J89 ( Supplementary  Figures S3A, S4C), indicates that drug BMS-202 did stabilize the PD-L1s dimerization system. To further understand the influence of drug insertion on the stability of the PD-L1 dimer, the RMSF of homo residues of 5J89 and IS PP in the 150 ns MD trajectories were calculated respectively. From Supplementary Figure S6, we can see that the insertion of BMS-202 can stabilize the entire dimer structure by reducing the motion flexibility of the residue conformations in IS PP .
The dynamic cross-correlation matrix (DCCM) (Hunenberger et al., 1995) was calculated to discuss the motion correlation in the trajectories of the 5J89 system and the IS PP system ( Figure 5). Each matrix element is the motion correlation coefficient between the two residues, and its value is between −1 and 1.
It can be clearly seen from Figure 5A that this matrix can be divided into four sub-blocks, among them two deep red sub-blocks arranged along the main diagonal show the strong-positivecorrelation inner-chain motion of the two chains in 5J89 respectively, which is a kind of integral motion within each chain. In contrast, there are no diagonal sub-blocks in Figure 5B, as the integral motion correlation of the chains in IS PP is weaker. The reason of difference may be due to that the drug insertion makes the inner-chain motion correlation of the residues stronger, and the random vibrations of the atoms in the chain change into the overall consistent motions. That is to say, the BMS-202 can stabilize the IS PP

Kinetics of Inter-State Transformation
Investigate the Kinetics of "BMS-202 Inserting Process" by ASMD In order to investigate the kinetic feasibility of "BMS-202 inserting process," adaptive steered molecular dynamics (ASMD) simulation (Ozer et al., 2010;Ozer et al., 2012) was performed to pull BMS-202 out slowly along the cavity (Supplementary Figure S7) in a near equilibrium state to simulate the reverse process of drug insertion. The cavity formed upon the dimerization interface by the native cleft of PD-L1 monomer is important for the drug stabilization (Zak et al., 2016;Guzik et al., 2017). On the basis of the Jarzynski's equality (Ozer et al., 2010;Ozer et al., 2012), the potential of mean forces (PMF) profile (Supplementary Figure S8) displayed the free energy changes of the system for the departure of the BMS-202 from the cavity.
Based on the PMF results of ASMD and the previous free energy results of MM-PBSA, we constructed a free energy surface profile (Supplementary Figure S9A) to describe the BMS-202 insertion process. It's worth noting that, the final PMF value of the steering process is 33.10kcal/mol, containing not only the stabilizing energy 12.82kcal/mol by the drug insertion in Final State: Holo PD-L1s Dimerization Mode in 5J89 section, but also the resistance energy by the residue side chains in the cavity. Thus, the energy barrier generated by the geometric resistance of the protein cavity is approximately 33.10 − 12.82 20.28kcal/mol. From the free energy surface profile (Supplementary Figure S9A), the free energy of the transition state in the "BMS-202 inserting process" is −10.87kcal/mol. That is to say, the "BMS-202 inserting process" is feasible in the kinetic aspect, as its transition state free energy is not higher than the transition state of three-body separation (0kcal/mol).
In addition, the following phenomena are observed in the ASMD trajectory: During the pull-out process, the biphenyl structure of the head of BMS-202 is significantly twisted, which shows that the angle between the two planes of the biphenyl rings is gradually reduced from the initial of 60°to about 0°degrees to adapt to the protein cavity. Moreover, the side-chains conformations of the "hydrophobic shrimp feet," "alanine clips" and "tyrosine sucking disks" of the protein interface are also twisted by the displacement of BMS-202. All of these phenomena reveal the key interactions between BMS-202 and PD-L1 dimers from another perspective. In combination with the previous analysis in BMS-202 Can Stabilize the IS PP section, it is reasonable to conclude that the insertion of BMS-202 produces a hydrophobic geometric complementary stabilization.
Compositing the transition state of "BMS-202 inserting process" (−10.87kcal/mol), the transition state of "three-body separation" (0kcal/mol) and the previous free energy results of all the stable PD-L1 dimerization modes by MM-PBSA, we constructed a integral free energy surface (Supplementary Figure S9B) to investigate the integral inter-state transformation kinetics in this multiple-states system.

Investigate the Kinetic Characteristics of Inter-State Transformation by GITR Theory
The inter-state transformation of three stable dimerization modes in the PD-L1 system can constitute a network (Figure 7), in which the existence of each dimerization mode in the system will affect the concentration and transformation rate of a specific dimer mode in the network. For this reason, the inter-state transformation kinetics is different from the traditional simple reaction kinetics problems. This kind of multiple-states self-assembly networks ( Figure 6) also exist in other systems (Hagan and Chandler, 2006;Rha et al., 2009;Pashuck and Stupp, 2010;Adamcik et al., 2011;Martinez-Avila et al., 2012;Grzybowski et al., 2017;Michaels et al., 2018;Dai et al., 2019). We attempt to provide a unified solution-generalized inter-state transformation rate (GITR) theory to this kind of problem.

Generalized Inter-State Transformation Rate Theory
For a complex self-assembly network containing any finitely many stable states (self-assembly modes) and any finitely many transition states (the Gibbs free energy maximum point on the transformation path), with arbitrary links between them, the transformation rates between self-assembly states in the network can be calculated by GITR theory. Our GITR theory has the following three basic hypotheses. To self-assembly system, they are rational and physical. 1) Single transition state hypothesis: It is reasonable because if a transformation between two stable states via two transition states (Gibbs free energy maximum points), we can find and insert a new stable state between the two transition states as a new addition. 2) The transition states are high energy and low concentration.
3) Inter-state transformation of the self-assembly system follows first-order kinetics. It is reasonable in many inter-state transformation problems of complex self-assembled systems (Odde et al., 1995;Cheng et al., 2011;Saric et al., 2016;Michaels et al., 2018). Because in many cases, the rate-FIGURE 6 | Each stable state S refers to the stable polymerization state formed by the aggregation of ν s self-assembly monomers. Each transition state T refers to the Gibbs free energy maximum point on the transformation path. There is no direct linking between any two stable states or any two transition states. For each transition state, we add a disk to distinguish it from the stable state, and the size of the disk is also measured by the relative Gibbs free energy. The R iti refers the transformation rate of stable state S i via transition state T i with the minus sign represents a decrease in concentration and the plus sign represents an increase in concentration. For the generality of the theory, the potential stable states and transition states are marked in gray. The C s which refers the concentration of stable state S is counted by polymerization cluster. The C 0 refers the total analytical concentration of selfassembly monomer in whole system. The G s and G ti are Gibbs free energies of the monomers in stable state S and transition state t i respectively. The Gibbs free energy of the free monomer state is set as the zero point of whole energy scale.
Frontiers in Chemistry | www.frontiersin.org November 2021 | Volume 9 | Article 783444 determining step of the transformation process from state S 1 to state S 2 of a self-assembly system depends on the energy obtaining of the monomer in state S 1 or even the conformation transformation to become the activated monomer in the S 1 -S 2 transition state (Cooper et al., 1983;Frieden, 1985;Kentsis and Borden, 2004;Jiang and Prentiss, 2014;Fu et al., 2015;Zierenberg et al., 2017). As a result, the rate-determining step is not intuitively the multi-monomer aggregation step.
Differential Rate Theory in the Network. From the classical transition state theory (CTST) (Eyring, 1935;Ross and Vlad, 1999;Fernandez-Ramos et al., 2006), the rate constant of elemental reaction is determined by the activation Gibbs free energy. As mentioned above, since the rate of inter-state transformation of many complex self-assembly systems is determined by the energy activation step of the monomer, the transformation rate from stable state S to transition state t i depends on the difference of monomer Gibbs free energy between the transition state t i and the stable state S. The transformation from the transition state t i to each directly linked stable state S is non-barrier process, so we propose the concept of linking fraction matrix L f . The sth row and t i th column element of L f represents the distribution ratio from the transition state t i to the directly linked stable state S. We will figure out the elements of linking fraction matrix L f later. Since the instantaneous transformation at a certain moment is an infinitesimal compared with C s , the transformation paths of a given stable state to any finitely many transition states will not affect each other. The transformation from stable state S to transition state t i is artificially selected to be negative and from transition state t i to stable state S to be positive. Thus, for any two directly linked states S and t i in the network, we get the corresponding path rate (15) where the C s which refers the concentration of stable state S is counted by polymerization cluster, ] s is the number of monomers in a cluster S, N A is Avogadro's constant, k B is Boltzmann constant, h is Planck's constant. All dummy indices summed over are marked with an apostrophe. Considering that in many cases, the activated monomer needs a dissociation process from the previous stable state or even a conformation transformation process, the rate difference between each path may be generated. (Cooper et al., 1983;Frieden, 1985;Kentsis and Borden, 2004;Jiang and Prentiss, 2014;Zierenberg et al., 2017) We introduce a dimensionless frequency factor A tis in the above formula to correct this process. Let k tis kBT h A tis and R N A k B , then Sum R st i over all the index t i that directly linked to the state S, then the following state rate R s represents the concentration change rate of stable state S: (17) Now let's calculate the path rate at equilibrium concentration. In order to obtain the path rate under equilibrium condition, we should firstly obtain the equilibrium concentration C eq s of the assembly state S. For the generality of the theory, we get C eq s by the probability P s that the monomers exist in the assembly state S at equilibrium.
Considering the grand canonical ensemble of self-assembly monomers under chemical potential μ, its grand partition function is where N is the number of monomers in the copy of the system with volume V, j is the jth quantum state of the system copy with volume V and number of monomers N, and E N,j is the energy of the jth quantum state of the system copy with number of monomers N.
Sum over the index j for a particular N, then the canonical partition function Q(N, V, T) of the system can be obtained. After substituting it in, the following equation can be obtained: Since the self-assembly monomer can be regarded as an indistinguishable delocalized particle, replace Q(N, V, T) by the canonical partition function of monomer q(V, T), then Sum over the index N, we get Calculate the grand canonical ensemble average for the number of monomers N: Since the volume V of all copies of the system in the grand canonical ensemble is constant, this indicates that the monomer concentration M is proportional to e μ k B T : We can write this relation in a more common form, for example, when self-assembly equilibrium is reached where M eq s and M°are respectively the concentration of monomer in self-assembly state S at equilibrium and the concentration of monomer in standard state. μ eq and μ°s are respectively the chemical potential of monomer at equilibrium and the chemical potential of monomer in the standard state of self-assembly state S. Therefore, the monomer concentration in each assembly state S at equilibrium is As a result, the probability P s that the monomers exist in the assembly state S at equilibrium is Note that, Z is constant for a given self-assembly system. Therefore, we get As a result, the path rate at equilibrium concentration is Since equilibrium has been reached and there is no net flow, the calculated path rate at equilibrium concentration should be 0: i.e., We can now see the physical meaning of linking fraction matrix L f : In a rough approximation, all the frequency factors A t i s can be equal. So all of the k t i s are equal to, let's say, k. Then the distribution ratio of the transition state t i to each stable state is equal to 1/n, where n is the number of stable states directly linked to the transition state t i . This is consistent well with the physical images. Because the transformation from the transition state t i to each directly linked stable state S is non-barrier process and all the k tis now is equal to k, and the k tis contains the information of dissociation frequency of reaction path, it means that all the dissociation frequencies of paths are equal, then the transformations of transition state t i to the direct linked stable states should be equiprobable, which are equal to 1/n. For the simple example shown in the figure above, its L f matrix should be where n s 6 and n t 3.
In the more precise case, L fst i k tis / s′ k t i s′ represents the transformation proportion of the transition state t i to the directly linked stable state S. It is determined by the probability of path selection from all the t i directly linked paths and the dissociation frequency along the reaction path together. For the simple example shown in the figure above, its L f matrix should accurately be s ′ 1 k t1s′ k t22 / 6 s ′ 1 k t2s′ 0 0 k t23 / 6 s ′ 1 k t2s′ 0 k t14 / 6 s ′ 1 k t1s′ k t24 / 6 s ′ 1 k t2s′ 0 0 k t25 / 6 s ′ 1 k t2s′ k t35 / 6 s ′ 1 k t3s′ k t16 / 6 s ′ 1 k t1s′ 0 k t36 / 6 s ′ 1 k t3s′ The single-direction interaction path rate at the equilibrium concentration can be calculated as follows: When we calculate the single-direction interaction path rate at equilibrium concentration in monomer form, it should be The above formula shows that the single-direction interaction path rate at equilibrium (in monomer form) is only related to the Frontiers in Chemistry | www.frontiersin.org November 2021 | Volume 9 | Article 783444 free energy of monomer in the transition state t i and the frequency factor A t i s of the corresponding path, but has no relation to the energy of stable state S. Again, this adheres to our intuition physics intuitive. The more stable the state with lower energy has a larger population, but it needs to climb a higher activation energy barrier. The transition state t i acts as the hinge of rate control in the network, and its monomer free energy G ti determines the probability of its occurrence in the system, so it determines the single-direction interaction path rate of the network.
Integral Time-Dependent Evolution of Concentrations in the Network. Since and dC s dt R s , write the above equation in vector and matrix form: In fact, both k t i s and the linking fraction matrix element L f st i record the topology information of the self-assembly network, because k tis and L fst i are 0 when the transition state t i and stable state S are not directly linked. After the above formula is simplified and combined, the following equation can be obtained: where ΔG represents the matrix k t i s e − G t i −Gs KT t i s nt×ns and δ sj is Kronecker notation.
Because we've already written the above formula in this dC dt AC form, C has to have a particular solution C ξe λt . After substituting it in, the following equation can be obtained: where E is the identity matrix. Since C ξe λt ≠ 0, ξ is the eigenvector of matrix A corresponding to the eigenvalue λ.
According to the superposition of the solutions, the general solution is where, z i is determined by the initial concentration condition, According to the integral results of the time-dependent evolution of the network, the concentration change rates of different stable states are determined by the same set of eigenvalues of the dynamic matrix A , and a variety of λ i can be regarded as a variety of time evolution scale factors. All the stable states in the network evolute in the same power law, the only difference is the front multiplicative factors which are the components of eigenvector ξ i . These front multiplicative factors can be regarded as the weight coefficients of each time scale in different stable states. This is also an interesting finding. In addition, if we let the matrix [k t i s ] evolve with the vector C and update by appropriate steps, the GITR theory would also be able to solve more complex non-linear kinetics problem.

Application of the GITR Theory to PD-L1 System
As a specific application of the GITR theory, we can now analyze the inter-state transformation kinetics in PD-L1 system. In Figure 7, the S 1 , S 2 and S 3 denote the initial state 1 (4Z18 apo form dimerization), the initial state 2 (IS PP apo form dimerization), the final state (5J89 holo form dimerization) respectively, T 1 and T 2 denote the three-body separation transition state and the "BMS-202 inserting process" transition state respectively. The relative free energy of the monomer in each state used in GITR theory is equal to the corresponding ΔG bind of state divided by the monomer number (according to Eq. 26), marked in black in Figure 7.
According to GITR theory proposed by us, under certain approximate conditions, there are unified path rate relations in the following form: Frontiers in Chemistry | www.frontiersin.org November 2021 | Volume 9 | Article 783444 R sti −k tis e − G t i −Gs RT × C s + s′ L fst i k tis′ e − G t i −G s′ RT ] s′ ] s C s′ (41) Then, the rates of the three transformation paths which are involved in the adjustment of PD-L1s from apo form dimerization (S 1 , S 2 ) to holo form dimerization (S 3 ) by BMS-202 are as follows where, for this system, the linking fraction matrix is L f ⎡ ⎣ k t11 / 3 s ′ 1 k t 1 s′ 0 k t 1 2 / 3 s ′ 1 k t 1 s′ k t 2 2 / 3 s ′ 2 k t 2 s′ k t 1 3 / 3 s ′ 1 k t1s′ k t 2 3 / 3 s ′ 2 k t2s′ ⎤ ⎦ ns×nt , and n s 3，n t 2.
For simplicity, we can further approximate that all the frequency factors A t i s are equal, hence all the k t i s kBT h A t i s are equal.
For this PD-L1s dimerization system, ] 1 2, ] 2 2, ] 3 2. Before the BMS-202 molecules are added, C 3 0, the concentrations of other states in the system should be the equilibrium distribution, hence C 1 and C 2 can be calculated by using the monomer distribution probability P s in GITR, where C 0 is the analytical concentration (total concentration) of PD-L1 monomer in the system: Then the rates of the three transformation paths above can be simplified as follows: where Z e − G 1 RT + e − G 2 RT . The minus sign represents a decrease in concentration, that is, a transformation to the final state S 3 .
When the values are substituted in, it can be found that the absolute value of R 2t 2 is the largest, that is, the corresponding path (S 2 to T 2 to S 3 ) is the dominant path for the dimerization mode transformation from apo form (pre-drug-adding) to holo form (post-drug-adding).

CONCLUSION
In this study, multi-scale computational simulation results and theoretical deduction systematically reveal the multiple dimerization mode transformation in the complex self-assembly system of PD-L1s with BMS-202.
In the thermodynamical aspect: Five states (two apo form dimerization modes, one holo form dimerization mode and two intermediate states) are identified as stable in the PD-L1 system. The binding site of the 4Z18 dimerization mode is mainly in the lower part (IgC domain) of PD-L1. The symmetry mismatch (between C 2v and σ v ) leads to the differences of binding free energies in 5J89 dimerization mode (between stripping mode 5J89 A and stripping mode 5J89 B ) and two intermediate states (between IS R and IS L ). A new stable apo form dimerization mode IS PP is found, which is more stable than the apo form dimerization 4Z18. The drug BMS-202 can stabilize the IS PP by a hydrophobic geometric complementary stabilizing effect which can reduce the conformational flexibility of PD-L1 interface residues and enhance the correlation of inner-chain residue motion.
In the kinetic aspect, we developed the GITR theory and then applied it to characterize the PD-L1 system. The results indicate that the BMS-202 can insert into and stabilize the apo dimerization mode IS PP (initial state 2) to form the final holo dimerization mode, and the "BMS-202 inserting process" is the dominant path of inter-state transformation in the PD-L1 system. Therefore, the drug BMS-202 acts more as a stabilizer than an inducer in the process of forming the holo form dimerization mode 5J89. Moreover, we give unified mathematical expression of the transformation rate of any path and state with an integral time-dependent evolutionary analytic solution in this type of complex self-assembly network which has arbitrary finite stable states, arbitrary finite transition states, arbitrary links between them, and arbitrary kinetic orders.
Our work would contribute to a better understanding of other complex multiple-states self-assembly network systems, both as an example and at the theoretical level.