Predicting the assembly/disassembly order of protein complexes via coarse-grained simulations

1.   Introduction

Protein complexes carry out the majority of catalytic, structural, and regulatory functions^[1]. Currently, cryo-electron microscopy (cryo-EM)^[2] and X-ray crystallography^[3] are the primary experimental techniques for obtaining high-resolution atomic structures of protein complexes. Additionally, computer techniques incorporating artificial intelligence enable the predictive modeling of protein complexes^[4]. With the continuous advancement of experimental and computational methods, an increasing number of structures of protein complexes have been determined, providing opportunities to investigate their functions at the molecular level.

However, in addition to static structures, it is essential to study the assembly of protein complexes^[5]. Many protein complexes have evolved to be assembled in a specific order^[6]. The stability and function of a protein complex rely on a multistep process of subunit assembly, which may provide valuable insights into the function and evolution of the protein complex. Numerous studies have demonstrated that improper protein assembly may result in negative biological effects^[7]. However, owing to the high complexity of the assembly pathway, the current understanding of the assembly mechanism is limited. According to Tompa and Rose’s study^[8], the assembly is hierarchical, and a protein complex must follow an ordered assembly pathway to achieve its biological function^[9]. Therefore, studying the assembly of protein complexes contributes to the reconstruction of protein complexes in vitro, the design of artificial protein complexes, and the development of new drugs related to assembly^[10].

The assembly of a protein complex can be studied by determining the assembly/disassembly order of its subunits. Stabilized intermediates of two or more subunits are detected via experimental techniques such as co-immunoprecipitation^[11], mass spectrometry (MS)^[12], and time-resolved cryo-electron microscopy^[13] or by creating deletion mutants to study the effects of deletions on assembly^[14]. MS is a powerful technique for identifying subcomplexes during the assembly/disassembly of a protein complex^[15]. A typical MS measurement generally starts with ionization of the protein complex^[16]. The complex is excessively charged, and signals from a set of subcomplexes are subsequently detected via MS. Subunits that bind each other weakly tend to dissociate from the complex^[17].

Computational modeling has the potential to predict the assembly/disassembly order of protein complexes. In early studies^{[18, 19]}, the order of subcomplex formation was predicted via a simple model based on the interface size, which was given by the number of contacted residues^[20]. With seven out of nine heteromeric complexes, the authors reported agreement between the interface sizes and assembly/disassembly orders. Path-LZerD^[21] uses knowledge-based potentials to evaluate binding preferences between subunits and can predict assembly paths even when the complex structure is unknown. The method involves constructing docked structures of subunits from single subunits via the Multi-LZerD method^[22], and iterative optimizations are performed via a genetic algorithm^[23]. The predictions of Path-LZerD are consistent with the experimental data for three complexes with up to four subunits, but for larger complexes with five or more subunits, the number of correct predictions decreased. The hybrid Monte Carlo/molecular dynamics simulation (hMC/MD) method^[24] uses a likelihood-based selection scheme to predict the disassembly order of a protein complex. The success rate of the method was greater than 0.9 for four tetrameric protein complexes.

When modeling large protein complexes, all-atom simulations are computationally expensive. Therefore, other methods, such as coarse-grained (CG) simulations, which are computationally efficient, have been proposed^[25–27]. Although a CG model does not describe protein interactions as accurately as an all-atom model does, the former can still capture the kinetic features of protein complexes^[26]. Furthermore, CG simulations can depict a biological process over a long timescale and have been employed for a variety of protein complexes^[27].

This work developed protocols for predicting the assembly/disassembly order of protein complexes via CG simulations. In one protocol, the disassembly process of a protein complex is simulated via CG simulations at different temperatures. This idea comes from MS experiments because standard electrospraying relies on the use of several factors, such as high temperatures. The temperature can be used to simulate the degree of ionization in MS and induce the dissociation of protein complexes. The disassembly temperature of each subunit is estimated through multitemperature CG simulations, after which the disassembly order of the complex is predicted. In the other protocol, iterative CG simulations of subcomplexes are carried out at room temperature to sequentially determine the most stable dimer, trimer, etc. Then, the assembly order of the protein complex is predicted by the stability of the subcomplexes in the CG simulations. Two protein complexes, the phosducin–Gtβγ complex and the Arp2/3 complex, were used to validate the protocols.

2.   Materials and methods

2.1   Protocol for predicting the disassembly order

To predict the disassembly order of multiple subunits in a protein complex, a protocol based on multitemperature CG simulations is used (Fig. 1). A high temperature disrupts interactions between protein subunits, which may lead to subunit disassembly. This process can be more easily observed by CG simulations than by all-atom simulations.

Figure  1.  Flowchart of multitemperature CG simulations for predicting the disassembly order of a protein complex.

DownLoad: Full-Size Img PowerPoint

The protocol consists of the following steps:

Step 1. The initial structure of the protein complex is prepared, such as by building missing regions.

Step 2. CG simulations are run at different temperatures. The temperatures are 300 K, 300 K + ΔT, 300 K + 2ΔT, …, until all the subunits are disassembled at a certain temperature.

Step 3. The CG trajectories are analyzed from low to high temperatures, and N_sub and ΔN_sub at each temperature are obtained. N_sub(T_i) is the number of disassembled subunits at the temperature T_i that can be obtained by analyzing the CG trajectory at T_i (details can be found in the next section). ΔN_sub(T_i) is the difference between N_sub at temperatures T_i and T_i−1.

Step 4. If ΔN_sub(T_i)=0, the temperature T_i is skipped. If ΔN_sub(T_i)=1, there is a disassembled subunit between T_i−1 and T_i, and we set the disassembled temperature of this subunit as T_i. If ΔN_sub(T_i)>1, there are two or more disassembled subunits between T_i−1 and T_i. In this case, ΔT_new needs to be set in the temperature range [T_i−1, T_i]. Similar to Step 2, a series of CG simulations are run at T_i−1+ΔT_new, T_i−1+2ΔT_new, …, until T_i−ΔT_new. These CG trajectories are analyzed as described in Step 3.

Step 5. The procedure was repeated until the disassembled temperature of each subunit was determined.

2.2   Determining the disassembly of a subunit via the q score

Native contacts are formed by pairs of amino acid residues that are physically close to each other in the native structure. A cutoff distance of the C_α atoms between two residues is set to determine if they form a native contact.

q score is the fraction of formed native contacts in a given conformation. In the native structure of a protein complex, the intersubunit native contacts are all counted. If two subunits have native contacts, the corresponding intersubunit q score is 1.0. During a CG simulation, the q scores may change. When the intersubunit q score between two subunits becomes 0, their native contacts are all broken. Therefore, q scores can be used to determine whether subunits disassemble.

2.3   Protocol for predicting the assembly order

To predict the assembly order of a protein complex, a protocol is developed by running iterative CG simulations of different subcomplexes, such as dimers and trimers (Fig. 2). The assembly order is determined according to the stability of the subcomplexes.

Figure  2.  Flowchart of iterative CG simulations for predicting the assembly order of a protein complex.

DownLoad: Full-Size Img PowerPoint

The protocol consists of the following steps:

Step 1. From the intact structure of the protein complex, all pairs of contacted subunits are extracted as the initial structures of dimers.

Step 2. CG simulations at room temperature (300 K) are run for all the dimers individually.

Step 3. The most stable dimer is found. For example, one can calculate the average root mean square deviation (RMSD) of each dimer during the CG simulation and pick the one with the smallest RMSD. The two subunits in this dimer are defined as core subunits.

Step 4. Starting from the most stable dimer, any other subunit in contact with the dimer is added to form a trimer. The initial structures of these trimers are taken from the protein complex.

Step 5. CG simulations are run for all the trimers individually.

Step 6. The most stable trimer is found.

Step 7. Steps 4 to 6 are repeated until an intact protein complex is formed. The assembly order is finally determined.

The initial structures of all the subcomplexes are taken from the intact protein complex, which means that we do not simulate the process of complex formation. If the CG simulations are started from separated monomers, the initial positions and orientations of the monomers affect their assembly order. Therefore, in this protocol, we predict the assembly order by simulating the stability of different subcomplexes.

Notably, this assembly protocol has several limitations. It does not consider the possibility that multiple dimers form first and then assemble from them, as well as other more complicated combinations. In principle, we could include all these possibilities, but there would be too many combinations of subcomplexes that require many CG simulations for a large protein complex.

2.4   The CG model

An atomic interaction-based CG model called AICG2+ was used for proteins^[28]. The energy function of the AICG2+ model is as follows:

R represents the 3n_aa Cartesian coordinates of the simulated protein, and R₀ represents those of the native structure. b_ibd is the ibd-th virtual bond length. r_ij is the distance between the i-th and j-th amino acids. $ {{V}}_{\mathrm{l}\mathrm{o}\mathrm{c}}^{\mathrm{f}\mathrm{l}\mathrm{p}} $ is a generic flexible local potential for virtual bond angles and dihedral angles constructed by analyzing loop structures in a protein structure database. The third term describes specific local interactions, and ε_loc,ij and W_ij in the Gaussian function represent the strength and width of the local interaction between residues i and j, respectively. ϕ_ij in the fourth term is the dihedral angle formed by the four consecutive residues i to i+3, and W_ϕ,ij is the width of this local interaction. The fifth term is the structure-based native contact potential, and a default value of 6.5 Å is set as the cutoff distance to determine native contacts^[28]. The last term represents the sum of interactions between non-native contact residue pairs.

2.5   CG simulations

All CG simulations were run via the CafeMol 3.2 package, which is a general-purpose CG simulation software for simulating proteins and their complexes^[29]. We used the AICG2+ model, excluded volume and electrostatic interactions for protein complexes. The electrostatic interactions were calculated via the Debye–Hückel equation with a dielectric constant of 78 F/m and an ion concentration of 0.15 M.

For each system, energy minimization was initially performed at 300 K, running 1000 steps with the steepest descent method, followed by 2000 steps with the conjugate gradient algorithm. A production simulation was subsequently run using Langevin dynamics at a certain temperature. The time step was 0.4 Cafe-time (1 Cafe-time ≈ 49 fs). In total, 10⁸ steps were run for each simulation, and a snapshot was saved every 5000 steps. The neighbor list was updated every 100 steps.

The parameter go_unit determines the weight of the Gō potential in the energy function, and its default value of 1.0 may make the complex rigid. Usually, different go_units from 0.5 to 1.0 are tested by running CG simulations at 300 K. Finally, a relatively small go_unit is chosen, with which the complex is not rigid in the CG simulation but essentially remains stable. When the unbinding/melting temperature of a protein is experimentally available^[30], it would be possible to optimize the go_unit by matching the experimental temperature in CG simulations. In such a case, the disassembly temperatures of the subunits obtained via CG simulations may correspond to the actual temperatures; otherwise, they are only used to judge the disassembly order of the subunits.

2.6   Protein complexes for testing the methods

We selected two asymmetric heteromeric protein complexes from the 3D Complex server^[20], the phosducin–Gtβγ complex (Fig. 3a) and the Arp2/3 complex (Fig. 3b), for testing the protocols.

Figure  3.  Structures of the two protein complexes used to test the protocols. (a) The phosducin–Gtβγ complex. (b) The inactive state of the Arp2/3 complex. The ATP/ADP molecules were excluded from the CG simulations.

DownLoad: Full-Size Img PowerPoint

The phosducin–Gtβγ complex. Phosducin binds tightly to Gtβγ of the heterotrimeric G protein transducin, preventing the reassociation of Gtβγ with Gtα–GDP and thereby inhibiting the G protein cycle^[31]. The crystal structure of the phosducin–Gtβγ complex (PDB 1A0R) from the bovine retina was solved with a resolution of 2.80 Å^[32]. The complex consists of three subunits, phosducin and the Gtβγ dimer (Fig. 3a). The β subunit forms a seven-bladed β propeller around which the γ subunit is wrapped, with the N-terminal helices of the β and γ subunits forming a coiled coil. Phosducin is composed of two domains: a helical domain at the N-terminus and a mixed αβ domain at the C-terminus. The former covers the top of the β propeller, whereas the latter covers one of its sides^[33]. According to experimental data, the Gtβγ dimer can be formed in the absence of phosducin; therefore, the order of assembly is known^[34].

The Arp2/3 complex. This complex is an important regulator of the cytoskeleton and a factor that enhances the mobility of cancer cells, showing a strong correlation with cancer development^[35]. The crystal structure of the inactive Arp2/3 complex (PDB 1K8K) from bos taurus was solved with a resolution of 2.0 Å^[36], and the nonterminal missing residues were built via homology modeling^[37]. The complex consists of seven subunits, including two actin-associated proteins (Arp2 and Arp3) and five other proteins (named ARPC1–ARPC5 in descending order of molecular mass)^[36]. Both Arp2 and Arp3 are structurally similar to actin and may form dimers with the ability to bind actin filaments to form branches^[38]. ARPC2 and ARPC4 hold the dimer together by antiparallel binding of the C-terminal long helix (Fig. 3b). The α/β domain of ARPC4 interacts with Arp2, ARPC1, and ARPC5, whereas the two similar α/β domains of ARPC2 interact with only Arp3^[39]. ARPC3 is only in contact with Arp3^[40].

3.   Results and discussion

3.1   Disassembly order of the phosducin–Gtβγ complex

First, we obtained a matrix of native contacts of the phosducin–Gtβγ complex (Table 1) via CafeMol^[29]. Gtβ and Gtγ have the greatest number of intersubunit contacts, which is consistent with the experimental data showing that the two subunits serve as dimers of the transducin.

Table  1.  Intersubunit native contacts in the atomic structure of the phosducin–Gtβγ complex.

Subunit	Gtβ	Gtγ	phosducin	Sum
Gtβ c(339)		a327	308	b635
Gtγ (65)			3	329
phosducin (188)				311
aNative contacts between two subunits; bthe total number of intersubunit native contacts of the subunit; cthe number in parentheses is the number of residues of the subunit.

 | Show Table

DownLoad: CSV

CG simulations were performed at different temperatures to simulate the disassembly process of the complex (Fig. 1), with go_units of 1.0 and 0.65. In the CG simulations at 300 K, the go_unit of 1.0 indicates that the conformations of the protein complex are rigid, whereas the go_unit of 0.65 allows significant conformational changes while still preserving the complex structure. The results indicate that the different values of go_unit may not affect the disassembly order of the protein complex; thus, the following presents only those results with a go_unit of 0.65.

For each subunit i, a weighted average of intersubunit q scores (denoted as <Q_inter>_i) was computed as follows:

where (qscore)_ij is the q score between subunits i and j, (Con_nat)_ij is the number of native contacts between the two subunits. The temporal evolution of <Q_inter> during the CG simulation at each temperature was investigated.

In the CG simulation at 367 K (Fig. 4), the <Q_inter> of phosducin (gray) first decreases to zero, followed by Gtβ (blue) and Gtγ (red). The dissociation times of the three subunits are very close to each other. The structure of Gtβ is largely unfolded when phosducin disassembles, after which the β and γ subunits dissociate very quickly. Therefore, the predicted disassembly order of the phosducin-Gtβγ complex is phosducin and the Gtβγ dimer.

Figure  4.  Time evolution of <Q_inter> for every subunit in the phosducin–Gtβγ complex during the CG simulation at 367 K. Several snapshots are shown.

DownLoad: Full-Size Img PowerPoint

The above procedure was repeated three times independently. Although the disassembly temperatures were different, the same disassembly order of the subunits was obtained. The results indicate that the predicted disassembly order is statistically reasonable.

3.2   Assembly order of the phosducin–Gtβγ complex

The assembly order of the phosducin–Gtβγ complex was investigated via iterative CG simulations of subcomplexes at 300 K (Fig. 2), and three independent simulations were run for each subcomplex. First, CG simulations were performed to find the dimer with the most core subunits. A comparison of the RMSD values of the three possible dimers revealed that the Gtβγ dimer is the most stable, with an average RMSD value of 2.0 Å (Fig. 5). In contrast, the phosducin-Gtγ dimer has a very large RMSD value because there are few contacts between the two subunits (Table 1), and they dissociate even in the CG simulation at room temperature. Therefore, the β and γ subunits assemble first to form the core dimer, followed by phosducin to form the trimer. Notably, our predicted assembly order is the inverse of the disassembly order predicted by CG simulations at 367 K (Fig. 4), which supports the notion that the assembly and disassembly processes are generally reversible in protein complexes^[19].

Figure  5.  Stability of dimers during the assembly of the phosducin–Gtβγ complex. The mean and standard deviation of the RMSD value of each dimer were calculated by taking the average of three independent CG simulations.

DownLoad: Full-Size Img PowerPoint

3.3   Disassembly order of the Arp2/3 complex

Table 2 lists the matrix of native contacts of the Arp2/3 complex. ARPC2 and ARPC4 have the largest number of intersubunit contacts, which is consistent with the experimental data showing that ARPC2 and ARPC4 serve as the core subunits of the complex^[40].

Table  2.  Intersubunit native contacts in the atomic structure of the Arp2/3 complex.

Subunit	ARPC4	ARPC2	Arp3	Arp2	ARPC1	ARPC5	ARPC3	Sum
ARPC4 c(167)		a244	72	102	188	107	0	b713
ARPC2 (282)			213	0	2	0	0	459
Arp3 (414)				79	0	0	88	452
Arp2 (378)					37	98	0	316
ARPC1 (368)						80	0	307
ARPC5 (143)							0	285
ARPC3 (173)								88
aNative contacts between two subunits; bthe total number of intersubunit native contacts of the subunit; cthe number in parentheses is the number of residues of the subunit.

 | Show Table

DownLoad: CSV

The go_unit was also set to 0.65, and CG simulations at different temperatures (Fig. 1) were performed to simulate the disassembly processes of the Arp2/3 complex.

At 315 K (Fig. 6a), the <Q_inter> values of Arp3 (orange), ARPC2 (cyan), and ARPC3 (purple) are essentially greater than 0.8. However, <Q_inter> of ARPC1 (green) decreases to zero in the middle of the CG simulation, followed by ARPC5 (yellow). That is, ARPC1 and ARPC5 disassemble sequentially at 315 K (Fig. 7). The <Q_inter> values of Arp2 (Fig. 6a, red) and ARPC4 (Fig. 6a, blue) decrease since they have lost contact with the disassembled ARPC1 and ARPC5. Arp2 is unfolded after the disassembly of ARPC1, although the former is still in contact with the core (Fig. 7). In the CG simulation at 316 K (Fig. 6b), Arp2 is disassembled since its <Q_inter> decreases to zero. The <Q_inter> values of Arp3 and ARPC4 decrease accordingly because they have lost contact with the disassembled Arp2. The <Q_inter> values of ARPC2 and ARPC3 remain high. The subcomplex ARPC2-ARPC4-Arp3-ARPC3 is well folded at this temperature (Fig. 7). At 350 K (Fig. 6c), ARPC3 disassembled. Notably, ARPC3 has only 88 native contacts with Arp3 (Table 2). However, the disassembly of ARPC3 requires a relatively high temperature, indicating that the binding of ARPC3 to Arp3 is strong. The <Q_inter> of Arp3 further decreases at 350 K because its contacts with ARPC3 have broken, and that of ARPC2 starts to decrease. The subcomplex ARPC2-ARPC4-Arp3 is largely unfolded at this temperature (Fig. 7). Finally, in the CG simulation at 352 K, Arp3 disassembles, followed by the ARPC2-ARPC4 core (Fig. 6d). The core dimer (Fig. 7) is the last to disassemble.

Figure  6.  Time evolution of <Q_inter> for every subunit in the Arp2/3 complex during the CG simulations at different temperatures. (a) 315 K, (b) 316 K, (c) 350 K, and (d) 352 K.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Snapshots during the disassembly process of the Arp2/3 complex at different temperatures.

DownLoad: Full-Size Img PowerPoint

According to the above analysis, the disassembly temperatures (order) of the peripheral subunits are 315 K (ARPC1 and ARPC5), 316 K (Arp2), 350 K (ARPC3), and 352 K (Arp3). If two or more subunits have the same disassembly temperature, we determine their disassembly order according to their <Q_inter> plots (Fig. 6). The procedure was independently repeated three times, and the same disassembly order of the subunits was obtained.

3.4   Assembly order of the Arp2/3 complex

The assembly order of the Arp2/3 complex was investigated via iterative CG simulations of subcomplexes (Fig. 2), with each subcomplex repeated three times independently. Information obtained from disassembly simulations or other sources can be used in the assembly protocol to determine the corresponding subcomplexes, thereby significantly reducing the number of CG simulations needed. Since ARPC2 and ARPC4 are known to be the core subunits, in the first iteration, only the CG simulation of the core dimer was conducted. The dimer is stable, with an average RMSD of 2.1 Å.

In the second iteration, each of the peripheral subunits except ARPC3 was added to the core dimer, and the corresponding trimer was used to run a CG simulation. ARPC3 is not considered in this iteration because it has no direct contact with the core dimer. ARPC2-ARPC4-Arp3 is the most stable trimer, with the smallest average RMSD of 2.0 Å (Fig. 8a). Therefore, Arp3 may assemble first to form the core dimer.

Figure  8.  Stability of subcomplexes during the assembly of the Arp2/3 complex. The mean and standard deviation of the RMSD value of each subcomplex are shown from (a) the second to (d) the fifth iteration. For each subcomplex, the RMSD was calculated by taking the average of three independent CG simulations.

DownLoad: Full-Size Img PowerPoint

In the third iteration, Arp2, ARPC1, ARPC3, and ARPC5 were added to the ARPC2-ARPC4-Arp3 trimer, and a CG simulation was subsequently conducted for each tetramer. As shown in Fig. 8b, the average RMSD of the subcomplex ARPC2-ARPC4-Arp3-ARPC3 is the smallest, which indicates that this tetramer has the highest stability. ARPC3 may assemble after Arp3.

We repeated the above procedure until a complete protein complex formed. In the fourth iteration, the most stable five-subunit subcomplex is identified as ARPC2-ARPC4-Arp3-ARPC3-Arp2 (Fig. 8c). Arp2 may assemble after ARPC3. Finally, in the fifth iteration, ARPC5 may assemble before ARPC1 because the former can form a relatively more stable six-subunit subcomplex than the latter (Fig. 8d). ARPC1 may be the last assembled subunit. Notably, although ARPC1 always has smaller average RMSD values than does ARPC5 in iterations 2–4 (Fig. 8a–c), the latter has a smaller RMSD than does the former after the addition of Arp2 to the subcomplex (Fig. 8d). This may be because ARPC5 can form more native contacts with Arp2 (Table 2), thus increasing the stability of the subcomplex.

Therefore, the predicted assembly order of the Arp2/3 complex is ARPC2, ARPC4, Arp3, ARPC3, Arp2, ARPC5 and ARPC1. This is the inverse of our predicted disassembly order via CG simulations at different temperatures.

The predicted assembly order of the Arp2/3 complex is in agreement with data obtained from several in vitro experimental studies. Reconstitution experiments indicated that ARPC2 and ARPC4 form a structural core of the complex^[40], which defines the center of rotation and translation of two subcomplexes. One subcomplex consists of ARPC2, Arp3, and ARPC3, and the other consists of ARPC4, Arp2, ARPC1, and ARPC5.

3.5   Comparison with methods based on the interface size

In this section, we compare our results with those obtained from a simple method based on the interface size^{[18, 19]}. Starting from an intact protein complex, its subunits are dissociated iteratively, and at each step, the subunit with the smallest interface size is removed. The interface size between two subunits is defined by the number of contacted residues, and the map of the interface size of a protein complex can be found at the 3D complex server^[20].

On the basis of the interface size map of the phosducin-Gtβγ complex, it is possible to predict the disassembly order, which is the phosducin and the Gtβγ dimer. This result is consistent with the disassembly order predicted by our method (Fig. 4).

According to the map of the interface size of the Arp2/3 complex, a disassembly order from the core subunit ARPC2-ARPC4 is predicted: ARPC3, ARPC1, ARPC5, Arp2, and Arp3^[20]. Compared with the disassembly order predicted by our method, the only difference is that of ARPC3. The interface size of ARPC3 is the smallest among the seven subunits in the Arp2/3 complex, so it is predicted to be the first dissociated subunit. However, in our CG simulations, the disassembly temperature of ARPC3 is 350 K, which is significantly higher than those of ARPC1 (315 K), ARPC5 (315 K), and Arp2 (316 K) (Figs. 6 and 7).

4.   Conclusions

Protein complexes play crucial roles in nearly all biological processes. The objective of this work is to examine the order in which subunits assemble and disassemble within a protein complex. We propose a protocol for predicting disassembly processes through CG simulations at multiple temperatures. This technique allows for a comprehensive examination of the disassembly order, comprising definable intermediate subcomplexes. Another protocol is developed that enables prediction of the assembly order via iterative CG simulations at room temperature. Despite the simplicity of the CG model, the results obtained for the phosducin-Gtβγ complex and the Arp2/3 complex are in fair agreement with some previous experimental data.

Notably, the assembly protocol does not consider the possibility that multiple dimers form first and then assemble from them, as well as other more complicated possibilities. In such cases, there would be too many combinations that require many CG simulations for all the possible subcomplexes. We may address this limitation by combining the information from the assembly and disassembly protocols. Since the disassembly simulations include all the possibilities, the two protocols can complement each other, as the assembly and disassembly orders may generally be reversible^[19].

In this preliminary work, only two protein complexes were used to test the protocols. There are certainly some issues to be addressed. We gathered information on more protein complexes and predicted assembly/disassembly orders via our methods and methods based on the interface size. Are the assembly and disassembly orders of protein complexes always reversible? Do multiple assembly/disassembly pathways exist? How many different assembly/disassembly orders can be predicted between these methods? Can these differences provide new insight into the assembly/disassembly of protein complexes? In the future, we also hope to expand the applicability of this method to additional biomolecular complexes, such as protein‒DNA and protein‒RNA complexes, by taking advantage of proper DNA/RNA CG models.

Preprint statement

Research presented in this article was posted on bioRxiv prior to publication in JUSTC. The corresponding preprint article can be found at https://doi.org/10.1101/2024.01.24.576999.