Models in this study were based on several PDB entries outlined in the SI section. Models were parameterized with the Amber ff99SB-ildn force field chosen based on our previous findings for enzymatic systems. Diethyl-phosphorylated (DEP) serine was parameterized as a single residue by combining force-field Ser parameters with ones derived for a triethyl-phosphate moiety. It was done with the parameterize tool by tailoring GAFF2 torsions towards QM scans on the wB97X-D level of theory with aug-cc-pVDZ basis set similarly to what was done before. Charges were assigned with the RESP protocol. Other ligands and Ser modifications used in the study (Figure ) were not used in unrestrained classical MD and thus did not require tailored parameters. For topology placeholders, they were parameterized with the GAFF2 force field with acpype with AM1-BCC charges. Each system was placed in a cubic box with periodic boundary conditions and solvated with tip3p-FB. All crystallographic water molecules were retained while all added were manually filtered in case of incorrect solvation. Na + and Cl -ions were added to neutralize the net charge and reach 0.15M ionic strength. Systems were minimized with 5000 steps of steepest descent. The equilibration phase consisted of seven steps and was performed in 5 replicates for each BChE system and in 3 replicates for others. First, an NVT run of 500 ps was performed while positionally restraining heavy atoms by 1000 kJ/mol/nm 2 . Velocity rescale thermostat was utilized for temperature coupling under 300K. Then, for five rounds of NPT equilibration of 500 ps each, restraint strength was gradually decreased as follows: 1000, 1000, 500, 200, 100 kJ/mol/nm 2 respectively. The NVT stage and the first NPT stage used a flexible water treatment (no LINKS constraints) for finer position adjustment. Stochastic cell-rescale barostat was used for pressure control. 2 fs time step was used in all steps. The final step was an unrestrained NPT run of 100 ns.