CHARMM Force Field for Protonated Polyethyleneimine Titus Adrian Beu, * Andrada-Elena Ailenei, and Alexandra Farcas ¸ We present a revised version of our previously published atom- istic Chemistry at Harvard Macromolecular Mechanics (CHARMM) force field for polyethyleneimine (PEI). It is based on new residue types (with symmetric C N C backbone), whose integer charges and bonded parameters are derived from ab initio calculations on an enlarged set of model poly- mers. The force field is validated by extensive molecular dynam- ics simulations on solvated PEI chains of various lengths and protonation patterns. The profiles of the gyration radius, end- to-end distance, and diffusion coefficient fine-tune our previous results, while the simulated diffusion coefficients excellently reproduce experimental findings. The developed CHARMM force field is suitable for realistic atomistic simulations of size/ protonation-dependent behavior of PEI chains, either individu- ally or composing polyplexes, but also provides reliable all- atom distributions for deriving coarse-grained force fields for PEI. © 2018 Wiley Periodicals, Inc. DOI:10.1002/jcc.25637 Introduction The design and practical development of effective gene car- riers, featuring high transfection efficiency, specificity, and bio- compatibility, are central to many modern gene delivery protocols. [1–5] Widely used as a nonviral gene vector, polyethy- leneimine (PEI: [CH 2 CH 2 NH] n ) occurs in linear or branched configurations. If protonated (with the NH groups partially replaced by NH + 2 groups), PEI shows a considerable buffering capacity which enables the condensation of DNA into polyplexes via electrostatic interactions between the proton- ated units and the negative phosphate groups of DNA. Due to inherent difficulties associated with developing realis- tic atomistic or coarse-grained force fields (FFs) for polycations, the number of theoretical/computational papers dealing in detail with solvated PEI chains or DNA-PEI polyplexes is rather limited. Even though some updates of the CHARMM FF signifi- cantly improved the treatment of DNA, [6,7] in order to conclu- sively investigate DNA-PEI condensation, a reliable FF for PEI is still needed. In one of the first systematic computational studies on the formation of DNA–polycation complexes, Ziebarth et al. [8] employed the Amber gaff FF [9] (notably not specifically parametrized for PEI) and, acknowledging the charge distribu- tion around the DNA helix to be the key issue to understanding DNA condensation, optimized the partial charges by the restrained electrostatic potential (RESP) method [10] based on ab initio data. The same atomistic FF was used in subsequent investigations on the protonation behavior of solvated linear PEI [11] and also in a recent study on PEI–DNA and PEI–siRNA complexes. [12] The molecular dynamics (MD) studies of Choudh- ury et al. [13] on the solvation dynamics of linear PEI essentially employed the same Amber FF, without notable improvements. Equally starting from the Amber FF, the partial charges of PEI were derived in the studies of Kondinskaia et al. [14] from ab initio calculations on four model trimers by the RESP method. The MD simulations on solvated DNA–PEI complexes of Sun et al. [15] adopted residues by analogy from the CHARMM27 FF, [16] and the torsional parameters, identified to be important, were improved by fits to ab initio data. As a first step in devel- oping a coarse-grained MARTINI FF for modeling the complexa- tion of RNA, Wei et al. [17] developed an atomistic FF for polyethylene-glycol-grafted linear PEI based on the CHARMM General Force Field (CGenFF) [18] using a “divide-and-conquer” strategy applied to small polymer building blocks. The dihedral parameters, in particular, were optimized relative to ab initio potential energy scans by using the Force Field Toolkit (ffTK). [19] Aiming for a more realistic modeling of the size- and protonation-dependent behavior of PEI, we recently published a new CHARMM FF for linear PEI. [20] As a major difference with respect to previous parametrizations, along with the partial atomic charges, we consistently adjusted the whole set of bonded parameters (for bonds, angles, and dihedrals), not only the dihedral contributions. The quality of the parametrization was enhanced by a more comprehensive body of basic ab initio data used in the optimization procedure (carried out by means of the ffTK application), namely stemming from two PEI model tetramers. Defining residues with C C N back- bone, we actually implemented a generic nonprotonated resi- due type and two fractionally charged residue types, the latter being employed in pairs to model the unitary protonation charge, which, according to the ab initio charge distributions, extends beyond the limits of single PEI monomers. Another T. Adrian Beu, A. Ailenei, A. Farcas¸ University Babes¸-Bolyai, Faculty of Physics, Departmentof Biomolecular Physics, 1 Mihail Kog alniceanu Street, Cluj-Napoca 400084, Romania E-mail: titus.beu@phys.ubbcluj.ro Contract Grant sponsor: Executive Unit for Financing Higher Education, Research, Development and Innovation (UEFISCDI); Contract Grant number: PN-III-P4-ID-PCE-2016-0474 © 2018 Wiley Periodicals, Inc. WWW.C-CHEM.ORG FULL PAPER Wiley Online Library Journal of Computational Chemistry 2018 1