Computational protocol: Molecular Dynamics Simulation Study of the Selectivity of a Silica Polymer for Ibuprofen

Similar protocols

Protocol publication

[…] The MD simulations were performed with the GROMACS 5.0.4 [] package applying the OPLS-AA [] force field, including the enhancements proposed by our group for sol-gel reagents in a recent publication []. GROMACS is an open source software package widely used to perform MD simulations to simulate a great variety of systems; it is one of the best programs to perform MD simulations due to its high speed and reliability—in particular, through the new GPU support, which is able to speed up a MD simulation up to 10 times []. All systems under study contained water, methanol, the anionic form of Ibuprofen (the template, IBU, A), the silica trimer SI3 (C), its anionic form SI− (D), and the dual cyclic silicate trimer corresponding to a hydrolysed and condensed species derived from the cationic dehydroimidazolium ORMOSIL (DHI+, B). All these structures are depicted in . In a previous work, we studied the DHI+ species with a complete –OH hydrolisation in the silica rings. In this case, we focused our investigation on testing the affinity of another ORMOSIL compound that is likely to be present in the mixture.Due to this, we present here a new DHI+ species with a replacement in a 1:2 ratios in the silica rings (C). This change was made in order to study the selectivity of the DHI+ molecule that is likely to be present in the real mixture. In previous works, we studied the affinity of other DHI species towards different templates. Using this kind of representation, we are able to study and simulate the affinity of two new DHI+ molecules for the template in order to evaluate which one has the better affinity for the template molecule. All the nonstandard parameters have been described in previous reports [,] except those used for the two new DHI+, which were parameterised and validated for the first time in this paper. With regard to the atomic point charges for the DHI+ and IBU species, these were calculated using GAUSSIAN 09 [] in an OPLS-AA compliant manner, meaning that the geometry was first optimised at the HF/6-31G* level, and partial charges were then computed from a single-point run, using the CHelpG scheme [] at the B3LYP/6-311++G(2d,2p) level of theory; information regarding these molecules is available under request. This approximation was chosen over the standard OPLS-AA force-field calculation (MP2/aug-cc-pVTZ//HF/6-31G*) due to a better stability of the DHI+ when using the 6-311++G(2d,2p) basis set. The studied model is specified in . The number of functional silicate (DHI+) units and structural silicate (SI3 plus SI−) units was determined from the experimental concentrations of DHI+ iodide and tetramethyl orthosilicate (TMOS). It was assumed that the precursors went through complete hydrolysis and fully condensed to DHI+ or SI3 (or its conjugated bases). On the other hand, the ratios of SI3 to SI− units and DHI+ were estimated from a species distribution analysis conducted at pH 9, based on the acidity constant of the silanol group, bearing in mind that pKa decreases by 1–2 units with high methanol contents.The initial state of the system was obtained using the PACKMOL package which inserts the respective number of units into the boxes at random positions []. Initial box dimensions were estimated considering the molecular weight and the density of each of the components of the mixture. After energy minimisation using the steepest-descent method implemented in the GROMACS package, a temperature annealing was performed in the NVT ensemble for 2 ns, attaining a temperature of 500 K, in order to ensure a proper mixing and gather three random independent initial configurations. These were, subsequently, used as starting configurations for the three independent MD equilibration runs needed to test the reproducibility of the simulations. Before the production stage, ~50 ns of simulation time in the NpT ensemble were taken to equilibrate the system and reach a stable configuration. Finally, production runs of 50 ns were performed in the NpT ensemble for data collection. Observable properties were sampled every 2 ps, from which total averages and standard deviations for each run were computed. The equations of motion were integrated using the Verlet leapfrog algorithm [], with a time step of 2 fs. Typically, the temperature (T) was kept fixed at 298 K by applying the velocity rescaling thermostat [], and, whenever necessary, the pressure (p) was held constant at 1 bar by using the Parrinello–Rahman scheme [,]. The time constant used for the Parrinello–Rahman coupling was set to 1 ps. Periodic boundary conditions were applied in all three Cartesian directions. For the water molecules, the Transferable Intermolecular Potential four-point model (TIP4P) [] was applied. The non-bonded electrostatic interactions were calculated using a sixth-order particle mesh Ewald (PME) method [] beyond a cutoff radius of 1.1 nm. The Lennard–Jones was calculated within a cutoff radius of 1.1 nm with the help of a neighbour list, updated every 10 time steps. A dielectric permittivity, εr, equal to 1.0 was used. Statistical and trajectory analysis of the simulations were carried out by resorting to the utilities included in GROMACS, while visualisations were made with Visual Molecular Dynamics (VMD) []. The analysis consisted essentially in the calculation of radial distribution functions (RDF), diffusion coefficients (D), and coordination numbers (NB), along with clustering analyses. The RDF between different types of molecules were calculated as (1)gAB(r)=〈ρB(r)〉〈ρB〉loc where 〈ρB(r)〉 refers to the average density of particle B at a distance r, around the particle A, and 〈ρB〉loc refers to the density of the particle B averaged over all spheres around particles A until a maximum radius (rmax), i.e., half of the box length. The RDF are additionally averaged on all particles of type A present in the system and averaged over the trajectory (simulation time). The g_rdf function included in the GROMACS package calculates the RDF in different ways. The normal method is around a (set of) particle(s), the other methods are around the centre of mass of a set of particles or to the closest particle in a set. Here, the RDF were calculated using both. The NB of a particle or atom B around another one A were calculated by integrating the radial distribution function between the centre of A and the first local minimum, rm: (2)NB=4πρB∫0rmgAB(r)r2dr where ρB refers to the density of species B (expressed in units of molecules per volume).The cluster analysis was performed using the g_cluster package included in the GROMACS software. This utility can cluster structures using several different methods. We determined structures from the trajectories of the runs using the single linkage, which adds a structure to a cluster when its distance to any element of the cluster is less than the cutoff. We performed the cluster analysis using cutoff values (i.e., the largest distance to be considered in a cluster) between 0.3 and 1.5 nm.As to the diffusion coefficients of the mixture components, these were calculated from the Einstein relation (mean-square displacement (MSD)): (3)D=16t〈|r→i(t)−r→i(0)|2〉 where r→i are the centre of mass positions of the molecules. The MSD is averaged over molecules, and, in order to improve the statistics, several restarts r(0) were used along the trajectory.Finally, the assessment of the validity of the potential considered for the simulations was demonstrated in our previous paper [] All the molecular species had been previously validated, save for IBU, which was parameterised for the first time here. The validation of this new molecule was performed just as in our previous work []. That is, a pure IBU system was considered and simulated, consisting of 100 molecules and 100 counter ions in a cubic box. The density of this template is 1.0 g/cm3, while the density calculated by MD was 1.084 g/cm3. These values appear to be acceptable according to the modifications and therefore, it allowed us to conclude that the parameterisation of the molecule is acceptable. […]

Pipeline specifications

Software tools GROMACS, VMD
Application Protein structure analysis
Diseases Genetic Diseases, Inborn
Chemicals Ibuprofen, Silicon Dioxide