Domenico Marson, Suzana Aulic, Maurizio Fermeglia, Erik Laurini and Sabrina Pricl
Reading time required
minutes
Speed reading only takes 22 minutes
š
Abstract
Nowadays, computer simulations have been established as a fundamental tool in the design and development of new dendrimer-based nanocarriers for drug and gene delivery. Moreover, the level of detail contained in the information that can be gathered by performing atomistic-scale simulations cannot be obtained with any other available experimental technique. In this chapter we describe the main computational toolbox that can be exploited in the different stages of novel dendritic nanocarrier productionāfrom the initial conception to the stage of biological intermolecular interactions.
In modern dendrimer-based nanocarrier design and development, computer simulations have been established as an indispensable tool. The details of the information provided by atomistic-scale simulations exceed the capabilities of current experimental techniques. In this chapter, we will introduce the core computational tools and methods available at various stages of dendrimer nanocarrier productionāfrom the initial concept design to the study of intermolecular interactions.
š
1 Introduction
Today computer simulations constitute an essential tool for the design and development of new dendrimer-based nanocarriers for both drug and gene delivery. Undoubtedly, this success is linked to the exponential increase in affordable computer power and to the optimization/development of highly scalable computer codes able to run in parallel even on consumer-grade graphical processing units (GPUs). Thus, for instance, under the magnifying lens of all-atom molecular dynamics (AA-MD) simulations researchers are now able to follow the exact time and space behavior of a given nanovector in solution at a molecular level, and to characterize its interactions with a specific target. Moreover, the same set of computational techniques allows performing in silico experiments dealing with, e.g., the interactions of dendrimers with biological membranes and proteins, and the eventual aggregation and clustering of dendrimers themselves. In essence, the analysis that can be performed with data obtained from AA-MD simulations is limited only by the creativity and coding proficiency of the investigator, and can be tailored to fit the specificity of the system under study.
With the rapid increase in computational power and the widespread application of efficient parallel computing codes, computer simulations have now become the core tool for designing and developing drug and gene delivery nanocarriers based on dendritic macromolecules. Through all-atom molecular dynamics (AA-MD) simulations, researchers can observe the dynamic behavior of nanocarriers in solution at a molecular scale and deeply analyze their interactions with target molecules. These computational techniques can also be used to simulate the interactions of dendritic molecules with biological membranes and proteins, and even explore their own aggregation behavior. The analytical potential of AA-MD simulations is virtually limitless, depending on the innovative thinking and programming skills of researchers, and can be flexibly adapted to different research systems.
Based on our own experience in the field of computer-assisted cationic dendrimer-based nanovectors for siRNA delivery, in this chapter we will describe the main AA-MD computational toolbox that can be exploited in the different stages of novel dendritic nanocarrier productionāfrom the initial conception to its structural and chemico-physical characterization to the stage of biological intermolecular interactions. Specifically, we will present a protocol for poly(amino amine) dendrimers (PAMAMs) as a proof-of-principle easily adaptable to other cationic dendrimer families. This protocol is subdivided in three main steps: (i) dendrimer and siRNA models building and optimization; (ii) AA-MD simulation of the dendrimer per se and in complex with its siRNA cargo, and (iii) simulation data analysis and correlation with experimental observations. We will further exploit the Notes sections to present some technical suggestions as well as some practical examples concerning the application of the computational techniques along with their biological relevance for the benefit of both computational and experimental scientists.
Based on our experience in utilizing cationic dendrimer-based vectors for siRNA delivery, this chapter will comprehensively introduce an AA-MD simulation toolbox covering the entire process from concept design to structural and physicochemical characterization, and to the study of interactions with biological molecules. Specifically, we will use poly(amino amine) dendrimers (PAMAM) as an example to demonstrate an easily extendable simulation process. This process includes three main steps: (1) building and optimizing models of dendrimers and siRNA; (2) conducting AA-MD simulations of the dendrimer itself and its complex with siRNA; (3) analyzing simulation data and comparing results with experimental data. Additionally, we will share technical tips and practical application cases in the Notes section to help researchers in both computational and experimental fields better understand the biological significance of these methods.
š
2 Materials
2.1 Forcefield and Software
AA-MD simulations try to mimic the behavior of individual atoms with the aid of a computer. Assuming all atoms as rigid spheres, and given the position in space of all atoms in a molecular system, the force experienced by each atom is computed applying a specific energy function, called a forcefield (FF). FFs can differ for the functional form and parameter sets used to calculate the potential energy of a given system of atoms. Newtonās laws are then used to integrate the FF over time, ultimately yielding the dynamic evolution (coordinates) of the molecular system under study.
While the underlying fundamental theories are highly similar, several different software (SW) are available to perform AA/CG-MD simulations. At the same time, a variety of FFs have been developedāfrom the most general ones to those custom-tailored to describe a particular systemāand not all of them are compatible with all current SW. In the specific case of AA-MD simulations of cationic-dendrimer nanovectors, the most widely employed FFs are the Chemistry at Harvard Macromolecular Mechanics (CHARMM) general forcefield, the General AMBER Forcefield (GAFF), the GROningen MOlecular Simulation (GROMOS) FF, and the Dreiding FF. Contextually, the principal simulation platforms adopted in the field are the Assisted Model Building with Energy Refinement (AMBER), the GROningen MAchine for Chemical Simulations (GROMACS), the Nanoscale Molecular Dynamics (NAMD), and the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS). Provided the availability of a given FF, the SW choice is mainly driven by the proficiency of the simulator with the different software and the accessible computational platform. For instance, the speed of AMBER is unparalleled when GPUs are used for computation, while LAMMPS has the steepest learning curve but offers high scalability on CPU-based supercomputers and an impressive variety of methods and procedures for simulation/analysis customization.
Although the fundamental theories are highly similar, there are many different software (SW) available for performing all-atom (AA) or coarse-grained (CG) molecular dynamics (MD) simulations, each with its own characteristics. Additionally, the existing forcefields (FF) range from general to specialized, but not all forcefields are suitable for all software. For AA-MD simulations of cationic dendrimer nanovectors, common forcefields include the CHARMM general forcefield, GAFF, GROMOS forcefield, and Dreiding forcefield. The main simulation platforms include AMBER, GROMACS, NAMD, and LAMMPS. The choice of software usually depends on the researcher’s proficiency with specific software and the computational platform being used. For example, AMBER performs exceptionally well with GPU-accelerated computing, while LAMMPS, although more difficult to learn, shows strong scalability on CPU-based supercomputers and offers extensive customization options for simulations and analyses.
The protocol reported in this chapter is based on the use of GAFF forcefield within the AMBER suite of programs. Our choice is based upon the following considerations: (i) the comprehensive availability of analysis tools within AMBER, (ii) the general agreement of the results obtained from simulations performed using GAFF FF with experimental data, (iii) its compatibility with AMBERās highly performing forcefields for biological macromolecules, and (iv) the unrivaled speed of AMBER simulation software on GPUs. In the used AMBER family of forcefields the potential energy function used assumes the form:
The protocol presented in this chapter is based on the use of the GAFF forcefield within the AMBER software suite. Our choice is based on the following considerations: (i) the comprehensive availability of analysis tools within AMBER, (ii) the general agreement of results obtained from simulations using GAFF FF with experimental data, (iii) its compatibility with AMBER’s efficient forcefields for biological macromolecules, and (iv) the unmatched speed of AMBER simulation software on GPUs. In the AMBER family of forcefields, the potential energy function used takes the form:
while the first three terms represent the interactions between two (bonded), three (angular), and four (dihedral) atoms connected via covalent bonds. The parameters kl,ākĪø, and Vn are force constants, reqand Īø~eq~ are equilibrium bond lengths and angles, respectively, and n and Ī³ are the multiplicity and the phase factor for the dihedral angles. The last term in this FF equation contains a standard Lennard Jones (LJ) term for the van der Waals interactions (in which rm, ij and Ļµ represent the minimum of the potential between two atoms i and j and the corresponding well depth, respectively), and a Coulombic energy term accounting for electrostatic interactions. In this term, qi and qj are the partial charges on the ith and jth atoms, and Ļµ~0~ is the vacuum permittivity. The rm, ij and Ļµij parameters for every couple of atoms in the system are computed from the self-interaction terms applying the LorentzāBerthelot rules: rm, ij = (rm, ii + rm, jj)/2 and
The first three terms represent the interactions between two (bonded), three (angular), and four (dihedral) atoms connected via covalent bonds. The parameters kl, kĪø, and Vn are force constants, req and Īø~eq~ are the equilibrium bond lengths and angles, respectively, and n and Ī³ are the multiplicity and phase factor for the dihedral angles. The last term in this forcefield equation includes a standard Lennard-Jones (LJ) term for van der Waals interactions (where rm, ij and Ļµ represent the minimum potential between two atoms i and j and the corresponding well depth, respectively), along with a Coulombic energy term that accounts for electrostatic interactions. In this term, qi and qj are the partial charges on the ith and jth atoms, and Ļµ~0~ is the vacuum permittivity. The parameters rm, ij and Ļµij for each pair of atoms in the system are computed from self-interaction terms using Lorentz-Berthelot rules: rm, ij = (rm, ii + rm, jj)/2 and
In summary, the first three terms describe interactions between two atoms (bonded), three atoms (angular), and four atoms (dihedral) that are connected by covalent bonds. The parameters kl, kĪø, and Vn represent force constants, while req and Īø~eq~ denote equilibrium bond lengths and angles, respectively. The last term includes a standard Lennard-Jones (LJ) term for van der Waals interactions (where rm, ij and Ļµ represent the minimum potential between two atoms i and j and the corresponding well depth), along with a Coulombic energy term that accounts for electrostatic interactions. In this term, qi and qj are the partial charges on the ith and jth atoms, respectively, and Ļµ~0~ is the vacuum permittivity. The parameters rm, ij and Ļµij for every pair of atoms in the system are calculated from self-interaction terms using Lorentz-Berthelot rules: rm, ij = (rm, ii + rm, jj)/2 and
2.2 Radius of Gyration and Asphericity
The radius of gyration Rg provides a quantitative characterization of the real size of a dendrimer molecule, and is defined as:
The radius of gyration Rg is an important quantitative indicator for characterizing the actual size of dendrimer molecules and is defined as:
where M and n are the total mass and number of atoms of the dendrimer, ci and mi are the position and mass of the ith atom, C is the center of mass of the dendrimer, and the angular brackets denote an ensemble average over the sampled configurations. The radius of gyration is strictly related to the spatial distribution of the dendrimer atoms, and ultimately to its size in solution. As such, Rg can be very useful to characterize the effect of some modification in the terminal units of a family of dendrimers of the same generation, or to study the effect of the protonation state on the structural configuration of a dendrimer. When computing Rg, it is also possible to estimate the relevant radius of gyration tensor Rg. The eigenvalues of Rg are its principal moments lx, ly, and lz (with lx ā¤ ly ā¤ lz), from the knowledge of which the dendrimer asphericity b can also be calculated as:
where M and n are the total mass and number of atoms of the dendrimer, ci and mi are the position and mass of the ith atom, C is the center of mass of the dendrimer, and the angular brackets denote an ensemble average over the sampled configurations. The radius of gyration is closely related to the spatial distribution of the dendrimer atoms and ultimately reflects its size in solution. Therefore, Rg is an important tool for assessing the effects of modifications in the terminal units of dendrimers of the same generation or for studying the impact of protonation states on the dendrimer’s structural configuration. When calculating Rg, the relevant radius of gyration tensor Rg can also be estimated. The eigenvalues of Rg are its principal moments lx, ly, and lz (with lx ā¤ ly ā¤ lz), from which the dendrimer asphericity b can also be calculated as:
As a useful descriptor of anisometry, b measures the deviation of a molecular geometry from a spherical form, in that for a perfectly spherical distribution of atoms b = 0, for prolate molecules b values are close to 0.25, whereas oblate molecules are characterized by values of b ā 1. As such, the evaluation of b can yield very insightful information about a dendrimerās shape modification when, e.g., it forms a complex with siRNA.
The asphericity b is a measure of how much a molecular geometry deviates from a spherical shape. For a perfectly spherical distribution of atoms, b = 0; for elongated molecules, b values are close to 0.25, while flattened molecules have values of b ā 1. Therefore, evaluating b provides valuable insights into the shape changes of dendrimers, for example, when they form complexes with siRNA.
2.3 Radial Monomer Density
While b and Rg are overall descriptors of a dendrimer shape, the radial monomer density Ļ(r) is a tool to study the distributions of monomers and/or other different molecular elements within a dendrimeric structure. Ļ(r) is defined as the number of atoms that are located within a spherical shell of radius r and thickness Īr from a reference center, usually taken as the dendrimerās center of mass. This analysis tool can convey many fine details information about the structure of the dendrimer under investigation including, by way of example, the extent of compactness of the dendrimer interior, the relative spatial distribution of the different dendrimer sub-generations within the entire dendrimer shell, and, last but not least, the distributions of the dendrimer terminal units with respect to the dendrimer core, i.e., the degree of dendrimer back-folding. Importantly, the distribution of other molecules like water, ions, counterions, and siRNA fragments with respect to dendrimer core can also be described by behavior of the corresponding Ļ(r).
While b and Rg provide overall descriptors of dendrimer shape, the radial monomer density Ļ(r) serves as a tool for studying the distributions of monomers and/or other molecular elements within a dendrimeric structure. Ļ(r) is defined as the number of atoms located within a spherical shell of radius r and thickness Īr from a reference center, typically the center of mass of the dendrimer. This analytical tool can provide detailed information about the structure of the dendrimer, such as the compactness of the dendrimer’s interior, the spatial distribution of different dendrimer sub-generations within the overall shell, and the distribution of terminal units relative to the core, indicating the extent of back-folding. Moreover, the distributions of other molecules like water, ions, counterions, and siRNA fragments relative to the dendrimer core can also be described through the behavior of the corresponding Ļ(r).
2.4 Solvent Accessible Surface Area and Interior Cavities
Water can play a key role both in the stabilization of the dendrimer structure and in its interactions with siRNA, mainly via the formation of an extensive network of hydrogen bonds. A measure of a dendrimerās surface exposed to the solvent is given by its accessible surface area (ASA). For any given molecule, ASA is typically calculated using the so-called rolling ball algorithm, which employs a sphere of a particular radius rp (e.g., a typical rp value is 1.4 Ć , which approximates the radius of a water molecule) to āprobeā the surface of the molecule. It is interesting to observe that, for a hypothetical spherical dendrimer without internal voids, the value of its ASA should increase with the square of the probe size rp
Water molecules play a crucial role in stabilizing the dendrimer structure and facilitating its interactions with siRNA, primarily through the formation of an extensive network of hydrogen bonds. The surface area of a dendrimer that is exposed to the solvent is measured by its accessible surface area (ASA). For any given molecule, ASA is typically calculated using the so-called rolling ball algorithm, which employs a sphere of a specific radius rp (e.g., a common rp value is 1.4 Ć , which approximates the radius of a water molecule) to āprobeā the surface of the molecule. Interestingly, for a hypothetical spherical dendrimer without internal voids, the value of its ASA should increase with the square of the probe size rp:
in which rd is the radius of the dendrimer. Accordingly, any deviation from linearity observed in a plot showing the square root of the calculated ASA as a function of rp is suggestive of the presence of internal voids in the corresponding dendrimer structure.
where rd is the radius of the dendrimer. Thus, any deviation from linearity observed in a plot of the square root of the calculated ASA against rp may indicate the presence of internal voids in the dendrimer structure.
2.5 MM/PBSA
A relatively quick and not particularly computationally intensive procedure to estimate the free energy of binding (ĪGbind) between a dendrimer and its cargo is the Molecular MechanicsāPoisson-Boltzmann (Generalized Born) Surface Area (MM-PB(GB)SA) methodology. The framework of this theory is summarized by the following equations:
A relatively quick and computationally inexpensive method to estimate the free energy of binding (ĪGbind) between a dendrimer and its cargo is the Molecular MechanicsāPoisson-Boltzmann (Generalized Born) Surface Area (MM-PB(GB)SA) methodology. The theoretical framework of this approach is summarized by the following equations:
where ĪEMM represents the gas-phase molecular mechanical energy change, composed by ĪEcovalent (the contribution of the covalent interactions energies: bonded, angles, torsions), ĪEvdW (the variation of the nonbonded van der Waals energy), and ĪEele (the change in the electrostatic calculated from the Coulomb potential). ĪGsolv represents the solvation free energy change, composed by its polar ĪGpolar and nonpolar ĪGnonpolar contributions. The sum of ĪEMM and ĪGsolv gives the enthalpic contribution to the free energy (ĪGbind), while TĪSbind is the corresponding variation in conformational entropy upon binding (where T is the system temperature in Kelvin). The polar term ĪGpolar is computed by replacing the solvent with a continuum medium with comparable dielectric constant (e.g., 78 for water) either using a Generalized Born (MM-GBSA) pairwise approximation or by directly solving the Poisson-Boltzmann eq. (MM-PBSA). The nonpolar contribution to solvation ĪGnonpolar, arising both from van der Waals interactions between the dendrimer/siRNA complex and the solvent and from the formation of a cavity in the solvent due to the presence of the solute itself, is usually calculated according to the following empirical expression for non-small molecules like dendrimers:
where ĪEMM represents the change in molecular mechanical energy in the gas phase, composed of ĪEcovalent (the contribution from covalent interaction energies: bonded, angles, torsions), ĪEvdW (the change in nonbonded van der Waals energy), and ĪEele (the change in electrostatic energy calculated from Coulomb potential). ĪGsolv represents the change in solvation free energy, which includes its polar ĪGpolar and nonpolar ĪGnonpolar contributions. The sum of ĪEMM and ĪGsolv provides the enthalpic contribution to free energy (ĪGbind), while TĪSbind reflects the change in conformational entropy upon binding (where T is the system temperature in Kelvin). The polar term ĪGpolar is calculated by substituting the solvent with a continuum medium of similar dielectric constant (e.g., 78 for water), either through a Generalized Born (MM-GBSA) pairwise approximation or by directly solving the Poisson-Boltzmann equation (MM-PBSA). The nonpolar contribution to solvation ĪGnonpolar, which arises from both van der Waals interactions between the dendrimer/siRNA complex and the solvent, as well as from the formation of a cavity in the solvent due to the presence of the solute itself, is typically calculated using the following empirical expression for larger molecules like dendrimers:
where ĪGdisp is the dispersion term due to the van der Waals interactions of the dendrimer/siRNA supramolecular assembly with the solvent, and can be computed via a solvent accessible surface integration. The rest of the equation accounts for the cavity formation via a combination of two empirical constants (Ī³ and Ī²) and a solvent accessible volume (SAV), i.e., the volume enclosed by the surface area of the solute. Finally, the entropic contribution to dendrimer/siRNA binding ĪSbind is usually estimated resorting to two most popular approaches: the normal mode analysis or the quasi harmonic approximation. As dendrimers are generally quite flexible molecules, a correct estimation of this last thermodynamic quantity may play a decisive role in the estimation of an ultimately realistic in silico ĪGbind value.
where ĪGdisp is the dispersion term due to van der Waals interactions between the dendrimer/siRNA supramolecular assembly and the solvent, calculated through solvent accessible surface integration. The remainder of the equation considers cavity formation through a combination of two empirical constants (Ī³ and Ī²) and a solvent accessible volume (SAV), which is the volume enclosed by the solute’s surface area. Finally, the entropic contribution to dendrimer/siRNA binding ĪSbind is generally estimated using two common methods: normal mode analysis or quasi-harmonic approximation. Given that dendrimers are typically quite flexible molecules, accurately estimating this last thermodynamic quantity may significantly influence the final realistic in silico ĪGbind value.