Our goal is to identify the stabilizing interactions that lead to formation of the C5b6 complex, the initial scaffold for the assembly of MAC. The crystal structure [32] shows three major sites of interaction between C5b and C6, named I, II, III (Figures 1A,B). Interactions between C6 and the thioester domain of C5b fall under Site I, while interactions of C6 with the macroglobulin (MG) ring of C5b fall under Sites II and III. The crystal structure also shows the presence of charged patches on the surfaces of the two proteins, C5b and C6 (Figures 1C,D), suggesting that charges may be contributing factors to binding. We performed an explicit solvent MD simulation, using a crystal structure as initial conformation, to optimize local geometries and chemistry and to delineate distinct conformational states visited throughout the simulation. Analysis of the MD trajectory, using the molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method, showed an overall favorable binding energy for the solvated complex, dominated by van der Waals interactions (Table 1). Given the observation of many charged patches on the surfaces of C5b and C6 (Figures 1C,D), but overall unfavorable polar contribution to binding (Table 1), we analyzed the frequency of occurrence of intermolecular pairwise polar interactions to obtain a closer look into the nature of polar contributions. Specific intermolecular salt bridges often stabilize protein complexes, an effect that has been termed ionic tethering, and so do intermolecular hydrogen bonds, even if the overall energetic contribution is dominated by van der Waals interactions of hydrophobic chemical groups [43].

We characterized the significance of intermolecular polar interactions by evaluating the number of salt bridges and hydrogen bonds that occur across the C5b6 binding interface at a frequency of at least 20% during the MD trajectory. Figures 2A,B shows occupancy (frequency of occurrence during the MD trajectory) maps for intermolecular salt bridges and hydrogen bonds, respectively. The C5b6 complex has 13 intermolecular salt bridges with distance cutoff of 5 Å, and 11 intermolecular hydrogen bonds, demonstrating varying levels of persistence throughout the trajectory.

The 13 intermolecular salt bridges (and protein/domain) in decreasing order of persistence (81–20%) are: Lys1117(C5b/C5d)-Glu2217(C6/CCP1), Lys884(C5b/CUB)-Asp2226(C6/CCP1), Glu646(C5b/MG Ring)-Arg1704(C6/LDLRa), Asp1076(C5b/C5d)-Arg2279(C6/CCP2), Lys1139(C5b/C5d)-Glu2187(C6/CCP1), Asp643(C5b/MG Ring)-Lys1693(C6/TSP2), Lys1139(C5b/C5d)-Asp2186(C6/CCP1), Asp648(C5b/MG Ring)-Arg1704(C6/LDLRa), Glu414(C5b/MG Ring)-Arg1734(C6/LDLRa), Glu149(C5b/MG Ring)-Arg2244(C6/CCP1), Asp648(C5b/MG Ring)-Lys1702(C6/LDLRa), Arg435(C5b/MG Ring)-Glu1716(C6/LDLRa), and Lys1133(C5b/C5d)-Glu2187(C6/CCP1) (Figure 2A). All intermolecular salt bridges (cutoff distance 5 Å) at all occupancy levels are listed in Supplementary Data Sheet 1, and all intermolecular Coulombic interactions (cutoff distance 8 Å, including salt bridges with cutoff distance of 5 Å) at all occupancy levels are listed in Supplementary Data Sheet 2.

Unlike the persistence of intermolecular salt bridges, most of the intermolecular hydrogen bonds formed between C5b and C6 appear less than half of the time in the trajectory. The 11 intermolecular hydrogen bonds (and protein/domain) in decreasing order of persistence (68–20%) are: Val1122(C5b/C5d)-Glu2188(C6/CCP1), Thr1105(C5b/C5d)-Thr2233(C6/CCP1) (backbone-backbone), Gly845(C5b/CUB)-Gln2241(C6/CCP1), Thr1105(C5b/C5d)-Thr2233(C6/CCP1) (side chain-side chain), Asn1077(C5b/C5d)-Arg2279(C6/CCP2), Asp1103(C5b/C5d)-Thr2233(C6/CCP1), Glu149(C5b/MG Ring)-Arg2244(C6/CCP1), Ser844(C5b/CUB)-Gly2239(C6/CCP1), Tyr41(C5b/MG Ring)-Asp1651(C6/TSP2), Ser96(C5b/MG Ring)-Pro1688(C6/TSP2), and His1106(C5b/C5d)-Arg2279(C6/CCP2) (Figure 2B). All intermolecular hydrogen bonds at all occupancy levels are listed in Supplementary Data Sheet 3.

Six salt bridges occur in Site III, followed by five in Site I, and two in Site II (Figures 1A, 2A). On the other hand, six hydrogen bonds occur in Site I, followed by three hydrogen bonds occurring in Site II and two hydrogen bonds occurring in site III (Figures 1A, 2B). Three bifurcated salt bridges are observed. These are Site I salt bridge Lys1139(C5b) with Asp2186(C6) and Glu2187(C6), Site III salt bridge Arg1704(C6) with Glu646(C5b) and Asp648(C5b), and Site III salt bridge Asp648(C5b) with Lys1702(C6) and Arg1704(C6) (Figures 1A, 2A). Two residues of C6 are hydrogen bonded with more than one partner, both in Site I. These are Thr2233(C6) side chain with Asp1103(C5b) side chain and Thr1105(C5b) side chain, and Arg2279(C6) side chain with Asn1077(C5b) backbone and His1106(C5b) side chain (Figures 1A, 2B). Finally, two residues of C5b and three of C6 participate in both salt bridge side chain-side chain interactions and hydrogen bonding side chain-backbone interactions. These are Glu149 and Asp648 of C5b, and Arg1704, Arg2244, and Arg2279 of C6. Among them are Arg1704(C6)-Asp648(C5b) salt bridge and hydrogen bond, involving Arg1704(C6) which is also part of a bifurcated salt bridge with Glu646(C5b) (Site III), Arg2244(C6)-Glu149(C5b) salt bridge and hydrogen bond (Site II), and Arg2279(C6) salt bridge with Asp1070(C6b), and hydrogen bond Asn1077(C5b) (Site I) (Figures 1A, 2).

Unfavorable charge-charge interactions within 5 Å of each other were not observed at occupancies greater than 10% (Supplementary Data Sheet 4). Thirteen unfavorable charge-charge interactions were observed within 8 Å of each other using the 20% occupancy threshold, as follows: four at Site I, three at Site II, and six at Site III. At Site I, the unfavorable charge-charge interactions are between Lys1139(C5b) and Lys1889(C6), Asp1140(C5b) and Glu2187(C6), Asp1167(C5b) and Glu2217(C6) and, Glu1181(C5b) and Asp2185(C6). At Site II, the unfavorable charge-charge interactions are between Arg177(C5b) and Arg2244(C6), Arg840(C5b) and Arg2252(C6) and, Lys884(C5b) and Arg2244(C6). At Site III, the unfavorable charge-charge interactions are between Asp43(C5b) and Asp1651(C6), Arg98(C5b) and Lys1690(C6), Arg435(C5b) and Arg1734(C6), Asp643(C5b) and Asp1706(C6), Glu646(C5b) and Glu1695(C6) and, Lys652(C5b) and Lys1702(C6). Of these, the only interactions with frequency above 50% are between Asp1167(C5b) and Glu2217(C6) at Site I, and Lys884(C5b) and Arg2244(C6) at Site II. All unfavorable Coulombic interactions with cutoff distances 5 Å and 8 Å at all occupancy levels are listed in Supplementary Data Sheets 4, 5, respectively. Overall, unfavorable charge-charge interactions are weak.

We further analyzed the MD trajectory to capture the dynamic changes in electrostatic contributions, by evaluating distinct conformational states. We identified six representative structures from the MD trajectory through principal component analysis (PCA) decomposition on phi and psi angles, followed by k-means clustering of the PCA components. The cluster centers (Figure 3) were extracted as representative structures. Structural representation of clustering conformational states is shown in Supplementary Figure 1.

We performed MM-PBSA calculations for the 20 most representative structures of each cluster to evaluate possible differences in the binding free energies among the six clusters. Although we observe differences among the clusters, the overall trends in their free energies are similar (Table 2).

We used the alanine scan method of our AESOP computational framework [44] to explore in detail the significance of the many pairwise charge-charge interactions in contributing to the stability of the C5b6 complex. We performed computational alanine scan by mutating every ionizable amino acid residue to alanine, one at a time, to generate families of single mutants for C5b and C6, and their C5b6 complex. Subsequently, we performed Poisson-Boltzmann electrostatic calculations to obtain electrostatic contributions to the free energies of binding for the binding reaction of each mutant. We present the results in reference to the parent structure, as differences between the calculated electrostatic free energies of binding of mutant structures minus the electrostatic free energy of the parent structure (see Methods). We performed the calculations for the representative structures of all six conformational clusters, and the results for those residues that participate in salt bridges (shown in Figure 2A) are presented in Figure 4. Loss of binding upon mutation is denoted by an electrostatic free energy of binding with a positive value, whereas gain of binding is denoted by a negative value. Loss of binding for the mutant indicates that the mutated residue is a contributor to binding of the parent protein by forming favorable charge-charge interactions. Likewise, gain of binding for the mutant indicates that the mutated residue opposes binding of the parent protein by contributing to unfavorable charge-charge interactions.

A few mutants show perturbations greater than thermal energy of 2k_BT at room temperature, most notably Lys1139A of C5b, followed by Lys884A, Asp1076A, Lys1117A, and Lys1133A (Figure 4A), and nearly every mutant of C6, most notably Arg1704, with the exception of Glu1716A (Figure 4B). These significant charge interactions depicted by the mutations of Figure 4 have been identified in the salt bridge occupancy maps of Figure 2A. AESOP data for all ionizable amino acid replacements within 8 Å of the C5b6 complex interface are listed in Supplementary Data Sheet 6. The most notable charge interactions, involving Lys1139(C5b) and Arg1704(C6) as depicted by the data of Figure 4, form strong bifurcated salt bridges in Sites I and II, respectively (see above). Another charge interaction, but of lower strength, involving Asp648(C5b) also forms a bifurcated salt bridge in Site III (see above).

The complement system contains thioester domains, referred herein as C3d and C4d, that can covalently attach to different surfaces and tag the cellular species with complement complexes and fragments. However, the C5d thioester-like domain of C5 is distinct from C3d and C4d because it does not have the ability to covalently tag a surface, or to exist as a standalone domain, but has gained the molecular capability to instigate MAC formation [32, 33]. Despite all three thioester/thioester-like domains, C3d, C4d, and C5d, having similar structures, C3d is a molecular hub for reactions that mediate inflammatory processes such as opsonophagocytosis and B-cell activation, whereas C5d participates in the binding of C6 and acts as the cornerstone for the MAC assembly. Thus, to elucidate on how C3d and C5d achieve their versatile functions, we performed a comparative analysis of the physicochemical properties of C3d and C5d, also including a comparative analysis of C3d and C4d for completion. Earlier studies have identified the residues of C3d that interact with FH [29] and CR2 [28]. FH is a potent regulator of C3b, and C3b's convertase complexes, that inhibits opsonophagocytosis on host cells. The CR2-C3d complex is formed as part of the B cell receptor-coreceptor complex and is a link between innate and adaptive immunity.

First, we performed a sequence alignment of all three domains, C3d, C4d, and C5d, to compare the conservation of their structurally and functionally important residues (Figure 5). The percent identity between C3d and C5d is 29.2%, between C3d and C4d is 35.7%, and between C4d and C5d is 28.8%. The contact sites with FH and CR2 are spread across the sequence of C3d. Of the total 310 C3d residues, 87 residues are involved in contacts with FH and CR2 combined (union of red and blue horizontal boxes in Figure 5). On the other hand, only nine C5d residues are involved in contacts with C6 (green horizontal boxes in Figure 5), and from those nine only one is a common site with C3d, that of Asp1104(C5d) and Gly1179(C3b). From the 87 residues that are involved in C3d-CR2/FH interactions, only 19 are homologous to C5d residues (intersection of red/blue horizontal boxes and identical C5d and C3d residues, denoted by vertical blue blocks, in Figure 5).

Lastly, of the nine residues of C5d that are involved in binding to C6, four are charged residues (Asp1076, Lys1117, Lys1133, and Lys1139) that interact by forming salt bridges with residues in the CCP1 and CCP2 domains of C6 (Site I, Figures 1, 2A). Contact map analysis showed that Lys1137 forms strong bifurcated salt bridges (Figure 2A). Also, the AESOP analysis showed that out of these four charged residues, alanine perturbation of Lys1139, produces the highest loss of binding mutation, hence making it one of the most important contact sites (Figure 4). In addition, Val1122 of C5d also makes a strong a hydrogen bond with Glu2188 in the CCP1 domain of C6, with a 68% frequency of contact from the total trajectory (Figure 2). Thus, from the nine residues of C5d that are involved in contacts with C6, Lys1139, and Val1122 play the most important role by participating in salt bridges and a hydrogen bond.

We also performed comparative binding/conservation analysis for C3d and C4d. There are 19 conserved C4d residues with C3d that are also found in the binding interfaces of C3d with FH and CR2 (union of red and blue horizontal boxes in Figure 5). The majority of the conserved C4d residues (17 out of 19) overlap with the FH contact sites of C3d (red-only horizontal boxes in Figure 5), with one of them overlapping with FH and CR2 contact site and two additional ones overlapping with CR2 only (blue-only horizontal boxes in Figure 5). The side chains of the two C4d/C3d conserved residues that overlap with CR2 contact sites are negatively charged. On the other hand, 11 of the 17 C4d/C3d conserved residues that overlap with FH contact sites have hydrophobic side chains, including the residue that overlaps with both FH and CR2 contact sites.

We also performed a comparative analysis of the electrostatic properties of C3d, C4d, and C5d. Earlier works have studied the properties of a highly negatively charged concave surface on C3d that is involved in recognition and binding of complement regulators [29] and receptors [28, 41] as well as bacterial inhibitors [30, 31]. Figure 6A shows the C3d concave surface with the negatively charged patch, surrounded by a neutral rim. C4d also has a negatively charged patch but spans mostly the cavity's periphery (Figure 6C). C5d on the other hand, contains a slightly shallower but more extended cavity, with a sparsely negatively charged patch on its center and distinct positively charged patches at the edges of the surface (Figure 6E). The size and charge density of the concave cavity in C3d may be the determining factor for binding to CR2, a property that is not shared by C4d and C5d. We also examined the electrostatic properties of the internal thioester bond face (Figures 6B,D,F), focusing on the small region that is responsible for covalently tagging pathogen surfaces after hydrolysis of the thioester bond. Both, C3d and C4d, contain positively charged patches, but the patch on C3d is more extended (Figures 6B,D). Unlike C3d and C4d, the positively charge patch of the thioester bond moiety is lost in C5d, which has slightly negatively charged character. In addition, this region forms a crevice, which is absent in C3d and C4d, to assist packing of three hydrophobic residues of Phe2172, Ile2174, and Met2175 (part of C6's FSIM sequence motif), as stated in [33]. C6 contacts C5d through an elongated linker that contains the predominantly hydrophobic FSIM sequence motif at one end, and a polar/negatively charged DDEE sequence motif, proposed here, toward the other end, before it loops out and returns to form additional contacts with C5d. The occupancy of intermolecular contacts of FSIM residues throughout the MD trajectory is as follows (>20% for 5 Å between heavy atoms): Phe2172-Gln898, 85.7%; Phe2172-Glu899, 62.6%; Phe2172-Asn902, 40.9%; Phe2172-Phe938, 39.2%; Ser2173-Gln898, 36.3%; Ile2174-Ala894, 25.0%; Ile2174-Gln898, 77.6%. The DDEE sequence motif comprises residues Asp2185, Asp2186, Glu2187, and Glu2188. Asp2186 and Glu2187 participate in the high-occupancy bifurcated salt bridge with C5b's Lys1139; Glu2187 also participates in bifurcated salt bridge with C5b's Lys1133 and Lys 1139; and Glu2188 participates in hydrogen bonding with C5b's Val1122, which is the highest occupancy hydrogen bond of the complex (Figure 2). Both, FSIM and DDEE side chains interact with part of Site I that is depicted by the ellipse of Figure 6F.

Formation of the MAC pore is propagated through the cleavage of complement C5 that forms C5b and C5a. The smaller fragment, C5a, is an anaphylatoxin that mediates inflammation through activation of immune cells. The larger fragment, C5b, interacts with C6 to form C5b6. The C5b6 complex interacts with C7 and C8 to form the C5b78 complex. Subsequently, multiple C9 molecules combine with C5b78 to form the transmembrane pore C5b6789₁₈ (MAC). We have recently developed quantitative models to study the dynamics of the alternative [24] and the combined alternative and classical pathways [25] of the complement system. These models are based on describing the biochemical reactions of the complement pathways using ordinary differential equations and are parameterized using experimental kinetic data. We used our latest model that encompasses the alternative and classical pathways to gain a systems level understanding of the terminal cascade of the complement system, starting at fluid phase C5 and ending in MAC pore formation on cell membranes [25]. Figure 7 shows the concentration-time profiles of C5, C5b, C6, C5b6, and MAC pores. As C5 is consumed (Figure 7A), MAC is produced (Figure 7B). C5 is minorly consumed to form C5b (Figure 7C), and so does C6 (Figure 7D) when it combines with C5b to form the C5b6 complex (Figure 7E). C5b and C5b6 show similar time responses, but with different concentration magnitudes. The time response of MAC pore formation shows a short lag phase, followed by an accelerated production phase. Unlike complement products C5b and C5b6 that are mostly consumed in ~40 min, the profile of MAC pores shows a continuous production after the 90th min mark.

Our results present the prompt activation and propagation kinetics of the complement system. In spite of the terminal pathway taking root later in the complement cascade (after the cleavage of C3), our results show MACs can be generated in less than 10 min. Furthermore, our results show despite the presence of complement regulators acting on the terminal cascade and numerous biochemical steps from C5 to C5b, C5b6, C5b67, C5b678, C5b6789, and the polymerization of C9 to form MACs (C5b6789₁₋₁₈), the rates of C5 or C6 consumption may be good predictors for the rate of MAC formation (Figure 7). This is highlighted where an accelerated rate of consumption for C5/C6 are present within the first 30 min, and this rate subsequently is reduced after the 30th min. Similarly, but inversely related to the concentration-time profiles of C5/C6, the MAC production exhibits an accelerated phase within the first 30 min that also subsequently diminishes after 30 min. However, it should be pointed out that the rate of production in MACs is still higher compared to the rate of C5/C6 consumption after the 30th min. In summary, kinetics of C5 or C6 consumption may be sufficient to study the kinetics of MAC formation despite the presence of numerous biochemical steps from C5/C6 to MACs.

The three-dimensional cocrystal structure of C5b6 with code 4A5W [32] was obtained from the Protein Data Bank (PDB). C5b is represented as chain A, whereas C6 is represented as chain B. All missing residues of C5b and C6 were added with MODELLER [66]. The crystal structure of C3d with code 3OED [35] was obtained from the PDB. The structures of C4d and C5d were extracted from the crystal structures of C4b and C5b6, with PDB codes 5JTW [36] and 4A5W, respectively. Missing residues for C4d were modeled using SWISS-MODEL [67]. Structural visualization and comparisons were performed using Chimera [68]. Molecular graphics were generated using Chimera. Protonation states of histidine were assigned using PDB2PQR [69]. The following transformations are needed to convert residue numbering used in 4A5W [32] to residue numbering used in our study: for C5b residues 19-765, subtract 18 from the 4A5W residue numbers; for C5b residues 766-1676, subtract 96 from the 4A5W residue numbers; for C6 residues 22-934, add 1559 to the 4A5W residue numbers.

Sequences of C3d, C4d, and C5d were extracted from the PDB entries 3OED [70], 5JTW [71], and 4A5W [32], respectively. Sequence alignments of C3d, C4d, and C5d were performed using Clustal Omega [72] and visualized with Jalview [73]. Identification of C3d residues involved in binding with FH and CR2 were acquired from MD simulation analysis performed in previous studies [28, 29]. Identification of C5b6 intermolecular interactions was performed as outlined in the MD trajectory analysis section, below.

Initial minimization of structure in the absence of water was performed using NAMD and the CHARMM36 force field [74, 75]. Subsequently, the structure was solvated in a cubic TIP3P water box leaving a minimum margin of 12 Å between any protein atom and the cube boundary. Sodium and chloride counterions were added to the system to achieve 150 mM ionic strength and neutralize protein charges. Adding water and counterions increased the total system size to 747,620 atoms. The solvated structure was energy minimized by undergoing 50,000 steps of conjugate gradient energy minimization before heating from 0 to 310 K with all protein atoms harmonically constrained to their positions after minimization. Next, five equilibration steps were performed in which the first five steps for a total time of 7 ns. During the four stages of equilibration, all protein atoms were constrained at a force constant of 10, 5, 2, and 1 kcal/mol/Å, respectively. The final equilibration step was concluded by only constraining backbone atoms with a force constant of 1 kcal/mol/Å. Following equilibration, AMBER16 [76, 77] was used for the production run for 100 ns with the following conditions: periodic boundary conditions, Langevin temperature control, a non-bonded interaction cutoff of 12 Å, with SHAKE algorithm used for constraining hydrogen bonds, and an integration time step of 2 fs.

Characterization and visualization of intermolecular interactions was performed using CPPTRAJ, pandas, and seaborn [78–80]. Analysis of buried solvent accessible surface area upon binding and visualization were performed with MDTraj [81] and matplotlib [82], respectively over all 2,000 frames in the trajectory. MSMBuilder [83] and MSMExplorer [84] were used for clustering, visualization and extracting representative structures from the trajectory for electrostatic analysis. Salt bridges between C5b and C6 residues was calculated with custom R scripts in conjunction with Bio3D package [85]. A distance cutoff of 5 Å was used. CPPTRAJ was used to analyze hydrogen bonds formed between C5b and C6 over the course of the trajectory. For hydrogen bonds, the default distance cutoff of 3 Å was used between acceptor to donor heavy atom, and an angle cutoff of 135°. To extract representative structures, PCA decomposition was performed on the phi and psi angles observed throughout the trajectory to reduce to four principal components using MSMBuilder. The MiniBatchKMeans method in MSMBuilder was utilized to cluster the four principal components to six distinct clusters and cluster centers were extracted as representative structures.

MM-PBSA calculations were performed using a thermodynamic cycle that decomposes the calculation of the free energy of binding into molecular mechanics (MM) force field calculations in a state of low dielectric coefficient (ε = 2) and solvation calculations by transferring the proteins from the low dielectric coefficient environment to a high dielectric environment (ε = 80). A frame interval of 4 was chosen for the MM-PBSA calculations, and hence a total of 500 frames from a total of 2,000 were processed. In our case, the MM calculations include van der Waals and electrostatic free energies, but not covalent geometry energy (bonds, angles, torsions) or entropic effects. We use the one-trajectory approximation, according to which we separate the structures of the components of the complex from the structure of the complex without additional minimization and without performing separate MD simulations of the complex components. Given that we used the one-trajectory approximation, covalent geometry energies and entropic effects are expected to cancel out in the binding scheme C5+C6→ΔGbind0C5b6.

The solvation free energy calculations include electrostatic contribution according to Poisson-Boltzmann electrostatic calculations, and nonpolar contribution (cavity solvation) described by an empirical term based on loss of solvent accessible surface area upon binding [86]. The following equations describe the MM-PBSA calculations. (1)ΔGbind0=GC5b60-(GC5b0+GC60)

Where the Binding Free Energy According to MM Force Field Calculations is Given by (2)ΔEbind,MM0=ΔEMM,vdW0+ΔEMM,electro0, and electrostatic contributions of individual components, C5b and C6, and complex, C5b6, to solvation free energy are given by Poisson-Boltzmann (PB) free energy differences (3)ΔGsolv,polar0=GPB,ε=800-GPB,ε=20.

The overall solvation contributions to binding, including polar and nonpolar effects, and are given by (4)ΔΔGsolv0=ΔGsolv,C5b60-(ΔGsolv,C5b0+ΔGsolv,C60)+ΔGnonpolar0, and the MM-PBSA free energy of binding is given by (5)ΔGbind,solv0=ΔEMM,vdW0+ΔEMM,electro0++ΔΔGsolv0

The Alanine scan method in the AESOP (Analysis of Electrostatic Structures Of Proteins) python package [44] was utilized to perform a computational alanine scan on ionizable residues at the C5b6 interface, and to evaluate their electrostatic contributions to binding. Alanine scans were performed for each of the six representative structures of the MD-derived conformational states.

AESOP utilizes the Adaptive Poisson-Boltzmann Solver (APBS) [87] to calculate grid-based electrostatic potentials, which are converted to electrostatic free energies. The program PDB2PQR is used to pre-assign charges and atomic radii for each atom according to the PARSE force field [88, 89], as well as to convert the PDB format to PQR format used by APBS. The selection of parameters for AESOP calculations has been described before [90].

Electrostatic free energies of binding were calculated according to a thermodynamic cycle that is similar to the one used for the MM-PBSA calculations, described above, except that binding at the reference state is evaluated using Coulomb's equation instead of the molecular mechanics method. The following equations are used, as described previously [44] and in recent applications of AESOP on C3d and its ligands [28–30]: (6)ΔGbind,solv0=ΔEbind,Coulomb0-ΔGsolv,C5b60-(ΔGsolv,C5b0+ΔGsolv,C60) where, (7)ΔGsolv0=ΔGelectro,ε=800-ΔGelectro,ε=200

Electrostatic free energies of binding for the family of alanine scan mutants are generated as deviations from the electrostatic free energy of binding of the parent protein, according to: (8)ΔGbind0=ΔGbind,solv,mutant0-ΔGbind,solv,parent0.

APBS is used to calculate electrostatic potentials for the solvation steps, and the program COULOMB (part the APBS suite) is used to calculate binding at the reference state. For the solvated state, dielectric coefficients of 78.54 and 20 were used for solvent and protein interior, respectively, while for the reference state a dielectric coefficient of 20 was used to resemble that of the protein interior [90]. The ionic strength of the solvated state corresponded to monovalent counterions of 150 mM concentration (physiological ionic strength), whereas the reference state had zero ionic strength. For the electrostatic analysis of the representative structures of C5b6 the number of grid points and mesh dimensions were set to 321 × 257 × 257 and 282 Å × 245 Å × 239 Å, respectively. For the electrostatic analysis of C3d, C4d, and C5d, the number of grid points and mesh dimensions were set to 84 × 98 × 96 and 129 Å × 129 Å × 97 Å, respectively.

The dynamics of the terminal cascade of complement system activation were modeled using a previously developed mathematical model that describes the biochemical reactions of the alternative and classical pathways [25]. The output of the model is reaction rates in the form of concentration-time profiles for all complement system proteins, enzymatic cleavage fragments, and association complexes. The model consists of a system of 290 ordinary differential equations (ODEs) and 142 kinetic parameters. Equations, initial concentrations, and kinetic parameters can be found in the Supplementary Information of [25]. The system of ODEs was solved using the ode15s solver of MATLAB (Mathworks, Natick, MA).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.