Minimal Coarse-Grained Model for Immunoglobulin G: Diffusion and Binding under Crowding

Immunoglobulin G (IgG) is the most common type of antibody found in blood and extracellular fluids and plays an essential role in our immune response. However, studies of the dynamics and reaction kinetics of IgG–antigen binding under physiological crowding conditions are scarce. Herein, we develop a coarse-grained model of IgG consisting of only six beads that we find minimal for a coarse representation of IgG’s shape and a decent reproduction of its flexibility and diffusion properties measured experimentally. Using this model in Brownian dynamics simulations, we find that macromolecular crowding affects only slightly the IgG’s flexibility, as described by the distribution of angles between the IgG’s arms and stem. Our simulations indicate that, contrary to expectations, crowders slow down the translational diffusion of an IgG less strongly than they do for a smaller Ficoll 70, which we relate to the IgG’s conformational size changes induced by crowding. We also find that crowders affect the binding kinetics by decreasing the rate of the first binding step and enhancing the second binding step.


S1 Brownian dynamics simulations
We performed Brownian dynamics (BD) 1,2 simulations using the propagation scheme due to Iniesta and García de la Torre 3 : where t is time, ∆t is propagation step, r is 3N -dimensional vector of positions, N is number of beads, k B is Boltzmann constant, T is temperature, D is 3N × 3N diffusion matrix, F is 3N -dimensional vector of forces, and X is 3N -dimensional vector of mutually-independent random normally-distributed variables with zero mean and unit variance.Square root of a matrix is defined by its Choleski factor B, i.e., To find the values of D and F at time t + ∆t/2, simpler propagation scheme due to Ermak and McCammon 1 is used: As the D matrix changes slowly with time, for the sake of computational efficiency we recomputed D and B every 100 steps instead of every step, following other papers [4][5][6][7] .

S1.1 Hydrodynamic interacions
We computed the position-dependent diffusion matrix D in a generalized Rotne-Prager-Yamakawa (RPY) approximation [8][9][10][11] : , where a M ij = max(a i , a j ) and a m ij = min(a i , a j ).Finally: for Here η is fluid viscosity, rij is normalized vector connecting i-th and j-th bead, r ij is distance between them, and a i is i-th bead's hydrodynamic radius.
To account for slow decay of D with distance (∝ r −1 ij ), we used Ewald summation scheme proposed by Smith et al. 12 .We set the parameter α controling the convergence of the Ewald summation to √ π (default value in BD BOX).The maximal magnitude of both real (m cutoff ) and reciprocal (n cutoff ) lattice vectors was 2.

S1.2 Repulsive interactions
For two large spheres composed of small particles of radii equal to σ interacting with the repulsive component of the Lennard-Jones potential (LJ12) with energy ε LJ , the total interaction potential is an integral of LJ12 potential over the volumes of the two spheres.Its polynomial expansion in r ij was derived by Henderson et al. 13 and for small separations can be approximated as follows 4 : Parameters for repulsive interactions are: ε LJ = 0.37 kcal mol −1 , σ = 0.15 nm.We included repulsive interactions between any pairs of particles, apart from the IgG beads belonging to the same macromolecule.
To avoid large forces causing numerical problems, for separations between the macromolecules' surfaces below r min = 0.3 nm, we kept the magnitude of force fixed and equal to F (r min) .The upper cutoff for the repulsive interactions was set to 25 nm.Following Ando and Skolnick 4 , we accounted for macromolecules' roughness by multiplying potential from Eq. (S5) by: setting h = 0.94 nm for Ficoll70, 1.21 nm for S, 0.30 nm for H, 0.73 nm for A 1 , A ′ 1 , A 2 , and

S1.3 Simulation details
We performed BD simulations with a customized version of the BD BOX software 14,15 .The customization concerned allowing for overlaps between bonded beads and performing diffusion matrix computation and Choleski decoposition every 100 steps, instead of every step, to decrease the computation time.
BD was performed in cubic boxes of 85 nm × 85 nm × 85 nm, and we applied periodic boundary conditions in all three directions.Time step ∆t was 0.5 ps and simulations lasted for at least 15 µs.Temperature was 298.15K and dynamic viscosity was η = 1.02 cP.
The compositions of mixtures studied in this work (i.e., molar fractions and number of particles) are shown in Table S1.

S1.4 Analysis
In every analysis procedure we discarded first 1 µs of the trajectory to ensure that the system is equilibrated.

S1.4.1 Translational diffusion
We obtained translational diffusion coefficients from Time-Averaged Mean Squared Displacement (TAMSD): where ∆t = 10 ns is a window length, r i is position vector of i-th bead, N traj is the number of trajectories, and N steps is the number of steps in a trajectory.The long-time diffusion coefficients were obtained by averaging D/D 0 between 3 and 5 µs.Uncertainty of the D l due to sampling error was estimated by dividing the simulations into 5 subsets and computing standard deviation of the mean, while treating the subsets as independent "measurements".

S1.4.2 Rotational diffusion
We obtained rotational diffusion coefficients from Orientation Autocorrelation (OA) 16,17 : where ℓ is the orientation vector and the window length ∆t is the same as in Eq. (S7).We defined orientation by a vector connecting S with A 2 , but we also verified that alternatives: S-H and A 1 -A 2 , lead to similar results.OA decays exponentially with time, and the decay rate is related to the rotational diffusion coefficient D ).The curves are calculated for two time intervals used for fitting, as indicated on the plot.D r is expressed in terms of the rotational diffusion coefficient in infinite dilution (D r,0 ).Rotational diffusion coefficient depends only weakly on ϕ occ , but its value is sensitive to the time interval where the fitting is applied.terminated when one of two events happens: 1. one of the arms of IgG gets very close to the antigen hot spot (r i < 4.94 nm, i ∈ {1, 2}) , 2. the hinge gets very far from the antigen hot spot (l > l esc = 25 or 35 nm).
Then we use the expression obtained by Northrup, Allison, and McCammon 18 to express the bimolecular diffusion-controlled reaction rate constant k D using the percentage of reactive trajectories β for a given l 0 and l esc .Assuming that all collisions are reactive, the expression reads 18 : where Ω is probability that molecule with l = l esc gets back to l 0 .Assuming that particles at long distances are well approximated by noninteracting (neither hydrodynamically nor directly) spheres, Ω = l 0 /l esc .
We performed BD-NAM simulations for ϕ occ = 0, ≈ 7.5, and ≈ 15%.Compositions of the systems are gathered in Table S2.We applied periodic boundary conditions in x and y directions, but in the z direction, we enclosed the system between two impenetrable walls.them, i.e.: For every system, we ran 800-1000 simulations.Uncertainties due to sampling errors were estimated by dividing the simulations into five subsets and computing the standard deviation of the mean.The NAM results obtained using Eq.(S10) are presented in Fig. S3.
The details of BD simulations are as follows.The time step ∆t was 0.5 ps.BD was propagated with the Ermak and McCammon scheme (Eq.(S3), but with ∆t/2 → ∆t).The temperature T was set to 293.15 K, and dynamic viscosity to η = 1.005 cP.The crowder particles interacted with each other and with IgG only with a hard-sphere potential.

S2.2 Second binding step
Subsequently, we used BD simulations to inspect the kinetics of the second step of IgGantigen binding.We placed additional antigen on the surface in distance of 6, 9, and 12 nm from the first one.Compositions of the systems are gathered in Table S3.We measured how much time it takes IgG bound with one arm to one antigen to bind to the second antigen with

r 16 :Figure S1 :
Figure S1: Macromolecular crowding and rotational diffusion of IgG.(left) Orientation autocorrelation function vs. time for three values of the occupied volume fraction ϕ occ measured for the SA 2 segment of IgG.(right) Rotational diffusion coefficient D r as a function of ϕ occ extracted from OA using Eq.(S9).The curves are calculated for two time intervals used for fitting, as indicated on the plot.D r is expressed in terms of the rotational diffusion coefficient in infinite dilution (D r,0 ).Rotational diffusion coefficient depends only weakly on ϕ occ , but its value is sensitive to the time interval where the fitting is applied.

Figure S5 :
Figure S5: Second step of antibody-antigen binding: cumulative distribution functions.

Table S1 :
Composition of studied mixtures

Table S2 :
Compositions of BD systems used to investigate the first IgG-antigen binding step.