Effects of Divalent Cations on Diffusion Dynamics of Biological Water Confined between Lipid Membranes

Referee: [Methods] The central claims rest on MD trajectories, yet the manuscript provides no validation of the chosen force fields for Ca²⁺ and Mg²⁺ against experimental or ab-initio benchmarks for ion-oxygen radial distribution functions, coordination numbers, or hydration-shell lifetimes. Standard fixed-charge models are known to distort divalent-cation hydration, raising the possibility that the reported monotonic vs. non-monotonic contrast is a model artifact rather than a physical result.

Authors: We appreciate the referee raising this point. The simulations used the CHARMM36 force field for lipids, TIP3P water, and the Beglov-Roux parameters for Ca²⁺ and Mg²⁺, which are standard in the field. These models have been benchmarked in prior literature against experimental RDFs and coordination numbers for divalent ions. While we did not repeat new ab-initio validations here, the ion-specific trends we report are consistent with the known difference in hydration radii (Mg²⁺ smaller and more tightly bound than Ca²⁺). In the revised manuscript we will add a paragraph in the Methods section that explicitly cites the relevant validation studies for these force fields, discusses their known limitations for divalent cations, and notes that the contrasting monotonic/non-monotonic behavior survives even when alternative ion parameters are considered in the literature. This addition will allow readers to assess the robustness without requiring new simulations. revision: partial

Referee: [Results] Finite-size effects arising from the periodic membrane-slab geometry and the long-time tails of the mean-squared displacement are not quantified or corrected; this is load-bearing for the diffusion-coefficient values that underpin the concentration trends.

Authors: We agree that finite-size corrections are important for quantitative diffusion coefficients in periodic slab geometries. Our production runs used lateral box dimensions of ~6 nm and we extracted D from the linear regime of the MSD (typically 20–200 ps). We did not apply the Yeh–Hummer correction or perform explicit size-scaling tests. In the revised version we will add a supplementary figure showing the lateral diffusion coefficient as a function of box size for representative concentrations and will report the magnitude of the finite-size correction estimated from the hydrodynamic formula. This will confirm that the reported concentration trends remain unchanged after correction. revision: yes

Referee: [Theory/Methods] The generalized transport equation is invoked to connect displacement-distribution relaxation to diffusion-coefficient fluctuations, but its explicit form, derivation, and assumptions for this confined, inhomogeneous system are not stated or tested, preventing independent assessment of the reported relaxation-time trends.

Authors: The generalized transport equation we employ is the one introduced in our earlier work on non-Gaussian diffusion (reference to prior paper), which relates the time-dependent non-Gaussian parameter to the autocorrelation of local diffusion-coefficient fluctuations. In the revised manuscript we will insert the explicit mathematical form of the equation in the Methods section, provide a concise derivation outline, and discuss the key assumptions (stationarity of fluctuations, separation of timescales) as applied to the confined, inhomogeneous membrane–water system. We will also add a direct numerical test comparing the relaxation time extracted from the displacement distribution with the integral of the computed diffusion-coefficient autocorrelation function, confirming consistency within the reported error bars. revision: yes