Finite-Size-Corrected Rotational Diffusion Coefficients of Membrane Proteins and Carbon Nanotubes from Molecular Dynamics Simulations

We investigate system-size effects on the rotational diffusion of membrane proteins and other membrane-embedded molecules in molecular dynamics simulations. We find that the rotational diffusion coefficient slows down relative to the infinite-system value by a factor of one minus the ratio of protein and box areas. This correction factor follows from the hydrodynamics of rotational flows under periodic boundary conditions and is rationalized in terms of Taylor–Couette flow. For membrane proteins like transporters, channels, or receptors in typical simulation setups, the protein-covered area tends to be relatively large, requiring a significant finite-size correction. Molecular dynamics simulations of the protein adenine nucleotide translocase (ANT1) and of a carbon nanotube porin in lipid membranes show that the hydrodynamic finite-size correction for rotational diffusion is accurate in standard-use cases. The dependence of the rotational diffusion on box size can be used to determine the membrane viscosity.


Details on Hydrodynamic Finite-Size Correction for Rotational Diffusion
In the main text, we showed that the hydrodynamic stream function ψ(x, y) for rotatory 2D flow in the membrane under periodic boundary conditions is approximated well by the 2D Green's function ϕ(x, y) of electrostatics in 2D. Here we give further details on the expansion of ϕ(x, y) in terms of harmonic functions with square symmetry. These expressions apply to the special case of square-shaped simulation boxes in the xy plane of the membrane, with L the width of the box. We write the Green's function as ϕ(x, y) = − ln r + πr 2 2L 2 + k≥4 a k p k (x/L, y/L) (S1) where r 2 = x 2 +y 2 . Here we defined the arbitrary constant in ϕ(x, y) so that lim r→0 [ϕ(x, y)+ ln(r)] = 0. The lowest-order harmonic functions in the unit square [−1/2, 1/2] 2 with non-zero coefficients are p 4 (x, y) = x 4 − 6x 2 y 2 + y 4 (S2) p 8 (x, y) = x 8 − 28x 6 y 2 + 70x 4 y 4 − 28x 2 y 6 + y 8 (S3) p 12 (x, y) = x 12 − 66x 10 y 2 + 495x 8 y 4 − 924x 6 y 6 + 495x 4 y 8 − 66x 2 y 10 + y 12 (S4) p 16 (x, y) = x 16 − 120x 14 y 2 + 1820x 12 y 4 − 8008x 10 y 6 +12870x 8 y 8 − 8008x 6 y 10 + 1820x 4 y 12 − 120x 2 y 14 + y 16 (S5) These polynomials are harmonic, ∇ 2 p k (x, y) = 0, and have square symmetry. The coefficients a k can be determined, e.g., by matching the k-th order derivative with respect to x of ϕ(x, y) + ln r at x = y = 0 to the analytic Lekner sum given as eq 10 by Grønbech-Jensen. 1 In this way, one obtains expressions for the coefficients in terms of rapidly converging sums that evaluate to: As explained in the main text, the lines ψ(x, y) = cϕ(x, y) define the stream lines of the hydrodynamic flow field, with c a constant defined by the boundary condition on the rotating cylinder. A possible concern is that ϕ(x, y) is not strictly constant on any circle with finite radius. Therefore, we determined the maximum difference of ϕ(x, y) from the mean value ϕ θ along circles of different radii using the Lekner sums given by Grønbech-Jensen. 1 For hydrodynamic radii R H = L/10, L/5, and L/3, we found maximum relative differences max θ |ϕ(cos θ, sin θ) − ϕ θ |/ ϕ θ of 0.0035 %, 0.075 %, and 0.76 %. We conclude that even for relatively large cylinder radii, the deviations are small.

Analysis of the Mean Squared Displacement
In the following, we explore if and when a regime of regular diffusion has been reached. We first use the dilute systems to estimate the linear regime from log-log plots and to compare the linear fits to a coupled-diffusion model. We then confirm our results for both dilute and constant-density simulations via additional fits at a later lag time interval. Linear Regime of the MSD. To assess when a linear regime of regular diffusion has been reached, we show the MSD curves in a log-log representation ( Figure S1). The initial spread in orientation resulting from fast molecular-scale dynamics leads to deviations from regular diffusion, MSD(t) = 2Dt, instead exhibiting apparent sub-diffusive behavior, MSD(t) ≈ c t α with α < 1 and c > 0. However, as shown in Figure S1 B, adding a simple offset, MSD(t) = 2Dt + a, accounts well for the dynamics after t ≈ 2 ns.
Model for Coupled Diffusion. In an effort to account for the MSD curves also in the initial non-linear phase (but outside the ballistic regime), we considered a model of harmonically coupled diffusion. Specifically, we assumed overdamped Langevin dynamics of the apparent rotation angle θ with a diffusion coefficient D 1 , and of a hidden angle γ with diffusion coefficient D 2 . One can think of γ as the orientation preferred by the slowly relaxing membrane environment, and of θ as the actual orientation. In this model, θ and γ are harmonically coupled by a potential k(θ − γ) 2 /2 with k the spring constant in units of k B T .
For this model, the apparent MSD of θ, averaged over the hidden dynamics of γ, is with a small added offset a accounting for the earliest dynamics. At short times, 0 < k(D 1 + the apparent friction coefficient acting on θ is the sum of the friction coefficients on θ and γ, resulting in an apparent diffusion coefficient D lim = 1/(1/D 1 + 1/D 2 ); at intermediate times, the dynamics transitions from free to coupled diffusion.
We fitted all MSD curves of the dilute systems over the entire time window from 0 to 7 ns to eq S10, giving us the individual diffusion coefficients D 1 and D 2 as well as the limiting diffusion coefficient D lim for large lag times. The resulting fits describe the MSD almost perfectly as indicated by the two example curves in Figure S1. The average offset is a = 2.4(6) × 10 −4 rad 2 and the average coupling constant is k = 292(55) rad −2 .
We compared the resulting D lim to the diffusion coefficient D from linear fits in the range from 3 ns to 6 ns. As shown in Figure S2, D lim and D agree well, with the average deviation being 0.13% and the largest deviation being 1.68%. This indicates that the chosen fitting range of the linear fits is already representative for the long-time behavior of the MSD. D lim from coupled model Figure S2: Comparison of the diffusion coefficients from linear fits in the range 3 ns to 6 ns to the long-time limiting values for diffusion coefficients from fits to the model of harmonically coupled diffusion, eq S10, in the range 0 ns to 7 ns. Results are for ANT1 rotational diffusion in simulations with a single ANT1 protein in the membrane at different box widths.
Interestingly, all three diffusion coefficients from the coupled model follow the size dependence eq 15 from our hydrodynamic theory ( Figure S3). This result implies that the  that the analysis at the earlier regime, from 3 to 6 ns, offers a good tradeoff between small systematic errors (because the asymptotic regime of regular diffusion has practically been reached) and small statistical errors (because the scatter in the MSD increases substantially at longer times). Most importantly, the finite-size effects on rotational diffusion are robust, independent of the fitting range, as can be seen by comparing Figure 3 of the main text and Figure S4.    Table 4: Simulation details of the constant-density simulations (from previous work 2 ) and