Molecular Dynamics Simulations of Complex Mixtures Aimed at the

Subscriber access provided by ONDOKUZ MAYIS UNIVERSITESI

Article

Molecular Dynamics Simulations of Complex Mixtures Aimed at the Preparation of Naproxen-imprinted Xerogels Riccardo Concu, Martin Perez, M. Natalia Dias Soeiro Cordeiro, and Manuel Azenha J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/ci5004575 • Publication Date (Web): 08 Nov 2014 Downloaded from http://pubs.acs.org on November 12, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

1

Molecular Dynamics Simulations of Complex Mixtures Aimed at the Preparation of Naproxen-imprinted

2

Xerogels

3 4

Riccardo Concu1*, Martin Perez1, M. Natália D. S. Cordeiro1*, Manuel Azenha2.

5 6

1

REQUIMTE/2CIQ-UP, Department of Chemistry and Biochemistry, Faculty of Sciences, University of

7

Porto, Rua do Campo Alegre, 687, 4169-007 Porto, Portugal.

8 9

ABSTRACT

10

The main objective of this study was to simulate for the first time a complex sol-gel system aimed at

11

preparing the S-Naproxen imprinted xerogel with an explicit representation of all the ionic species at pH 9.

12

To this purpose, a series of molecular dynamics (MD) simulations of different mixtures, including species

13

never studied before using the OPLS-AA force field, were prepared. A new parameterization for these

14

species was developed and found acceptable. Three different systems were simulated representing two types

15

of pregelification models: the first one representing the initial mixture after complete hydrolysis and

16

condensation to cyclic trimers (model A), the second one corresponded to the same mixtures after

17

evaporation process (model B) and the last one (model C) a simpler initial mixture without an explicit

18

representation of all imprinting-mixture constituents. The comparison of systems A and C served mainly to

19

the purpose of evaluating if an explicit representation of all the components (model A) was needed, or if a

20

less computationally demanding system, where the alkaline forms of the silicate species were ignored (model

21

C), was sufficient. The results confirmed our hypothesis that an explicit representation of all imprinting-

22

mixture constituents is essential to study the molecular imprinting process, because a poor representation of

23

the ionic species present in the mixture may lead to erroneous conclusions or lost information. In general, the

24

radial distribution function (RDF) analysis and interaction energies demonstrated a high affinity of the

25

template molecule, 2-(6-methoxynaphthalen-2-yl) propanoic acid (S-Naproxen, NAP-), for the gel backbone,

26

especially targeting the units containing the dihydroimidazolium moiety used as functional group. Model B,

27

representing a nearly gelled sol where the density of silicates and solvent polarity were much higher

28

relatively to the other models, allowed for much fast simulations. That gave us the chance to observe the

29

templating effect through the comparative analysis and observation of the trajectories of simulation with the

30

template vs. non-template containing mixtures. Overall, a strong coherence between the imprinting-relevant

31

interactions, aggregation or the silicate network texturing effects taken out of the simulations and the

32

experimentally high imprinting performance and porosity features of the corresponding gels was achieved.

33 34

1. INTRODUCTION

35

In the past few years the sol-gel polycondensation technique has been increasingly employed with great

36

success as an alternative approach to the preparation of molecularly imprinted materials (MIM)1, 2, 3. Using

37

this approach it is possible to prepare synthetic host systems bearing improved molecular selectivity4. With

38

respect to classical techniques, MIMs present several advantages such as physical robustness, long shelf life,

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

39

simple preparation, great selectivity, etc. For these reasons, MIMs are being studied and used in very

40

different fields such as solid phase extraction, enantiomer separations, drug delivery, drug discovery and

41

ligand binding assays. These materials are often based on a silica backbone and on inorganic-organic hybrid

42

materials prepared with organically-modified trialkoxysilanes (ORMOSIL). The most widely used

43

precursors for preparing sol-gel materials have been the silicon alkoxides, such as tetramethoxysilane

44

(TMOS) or tetraethoxysilane (TEOS). Commonly, the imprinting protocols follow a co-polycondensation

45

route between one of these precursors and an ORMOSIL, R-Si(OR’)3 or (OR’3)Si-R-Si(OR”3), where R

46

represents an organic motif with affinity (usually via non-covalent interactions) for the template of

47

imprinting. Upon hydrolysis and poly-condensation of the precursor and ORMOSIL, a hybrid organic-

48

inorganic network is developed, where the [template-R motif] complexes are trapped upon gelification.

49

Finally, the removal of the template leaves a vacant recognition R-site, with the ability to recognize the

50

shape, the size and the functionality of the template (or structurally related compounds). In Figure 1 we show

51

a general scheme of the MIM process under study in the present work. It focuses on the dehydroimidazolium

52

motif, a cationic ORMOSIL (Figure 2A), aimed at sol-gel molecular imprinting of the carboxylate form of

53

the non-steroidal anti-inflammatory drug 2-(6-methoxynaphthalen-2-yl) propanoic acid (S-Naproxen, NAP-,

54

Figure 2E). We have shown5 that highly selective sol-gel naproxen imprints may be prepared by the co-

55

polycondensation between hydrolyzed TMOS and 1-(triethoxysilylpropyl)-3-(trimethoxysilylpropyl)-4,5-

56

dihydroimidazolium iodide(AO-DHI+ iodide) in a water and methanol mixed solvent, using NAP− as

57

template. During those studies, several experimental observations raised questions whose answers were not

58

straightforward. We could only appeal to our chemical intuition on the molecular scale phenomena when it

59

came to answer questions such as why the end-capped imprints were not so effective or why no

60

enantioselectivity could be observed. We believe that these, and many other macroscopic observations

61

related to the field of sol-gel molecular imprinting, may be scrutinized via the application of computational

62

chemistry techniques. In fact, we recently demonstrated, for the first time, the usefulness of molecular

63

dynamics (MD) simulations of pre-gelification sol-gel mixtures6,

64

aspects of molecular interactions, both related to the imprinting process or network structuring. Those MD

65

simulations did not model the complex reactive mechanisms involved in the sol-gel process, a too complex

66

task considering the current computational resources. Non-reactive mixtures, representing certain stages of

67

the sol-gel process, were thus considered as simplified approaches to the modeled systems. The mixtures

68

contained water, methanol, damascenone (the imprinting molecule), Si3O3(OH)6 (SI3, Figure 2C) and, some

69

of them, Si3O3(OH)5C3H6NHC6H5(SIPA). SI3 and SIPA are 3 silicate rings used as representatives for the

70

initial, small condensation products obtained from the polycondensation of 3 TMOS units (SI3) and the co-

71

polycondensation of TMOS with the ORMOSIL (3-propylaminophenyl)-trimethoxysilane at a 6.7:1 molar

72

ratio. We were conscious of such oversimplification of the real sol-gel mixture, and meant mainly to join

73

together (respecting the experimental proportions) species regarded as providing a good resemblance with

74

the final gel backbone (namely the same surface groups), the template and the solvent medium. All the

75

components were neutral, since the simulations respected to sol-gel synthesis were carried out at low pH

76

values (below silica’s isoelectric point, ~2). The simulations proved to be a useful tool for studying the

7

forward-gaining insight into the fine


Page 2 of 51

Page 3 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


77

solvation environment of the template and the silica clusters. The collected information about the competing

78

interactions in different simulated systems, corresponding to different experimental conditions, was in

79

excellent agreement with experimental observations of the performance of the prepared damascenone-MIMs.

80

Later, a similar study allowed us to ascertain the possible effects of adding polyethylene glycol to the sol-gel

81

mixtures used for the preparation of the damascenone-imprints3.

82

Herein we present a study aimed at extending the applications of MD simulations for the investigation of the

83

fine, atomistic level aspects related with the preparation of naproxen-imprinted hybrid silicas, as referred

84

above. In this case, better imprinting performances were experimentally observed when NAP−, the sodium

85

carboxylate derivative of naproxen, was employed as template, and the sol-gel process took place at pH ~9.

86

To be able to cope with such alkaline pH conditions, the MD simulations needed to incorporate a series of

87

charged species, in agreement with the distribution predicted for the several species at that pH, according to

88

the acid-base equilibria involved. Moreover, the organic motif, targeting the imprinting of NAP− is a

89

permanent cation of the (OR’3)Si-R+-Si(OR”3) type, a more complex structure than used before. Additionally,

90

counter anions and cations will also be needed in the system. Comparing to the previous MD simulations,

91

these are of higher complexity and justify an effort to demonstrate the feasibility of MD under such “high

92

charge density” conditions. That constitutes the main goal of the present research, which may be regarded as

93

an intermediate, but necessary, stage before the envisaged application to more advanced simulations, such as

94

the sol-gel/oil interfacial studies in the context of the emulsion preparation of imprinted silica micro- or

95

nano-spheres.

96 97 98

2. MATERIALS AND METHODS

99

2.1. Computational Details. The MD simulations were performed with GROMACS 4.5.58 package

100

applying the OPLS-AA9 force field, including the enhancements proposed by Price et al17 for sol-gel

101

reagents. GROMACS is an open source software package widely used to perform MD simulations to

102

simulate a great variety of systems; it has a good reputation concerning speed and reliability, in particular

103

through the new GPU support which is able to speed up a MD simulation up to 10 times. All systems under

104

study contained water, methanol, the carboxylate form of S-Naproxen (the template, NAP−), the dual cyclic

105

silicate trimer (Figure 2A) corresponding to a hydrolyzed and condensed species derived from the cationic

106

dehydroimidazolium ORMOSIL (DHI+, [Si3O3(OH)5)2-(C3H6)]2-C3H5N2+, and the silica trimer SI3.

107

Depending on the case, some contained also DHI+/− (the zwitterionic form of DHI+, Figure 2B), obtained

108

from DHI+ by the loss of a proton, and also the anionic form of the silica trimer (SI3−, Si3O3(OH)5O−, Figure

109

2D). All the nonstandard parameters were described in a previous report10, except those used for DHI+,

110

DHI+/− and SI3−, which were parameterized for the first time in this work. With regard to the atomic point

111

charges for the DHI+, DHI+/− and SI3− species, they were calculated using GAUSSIAN 0911 in an OPLS-AA

112

compliant manner; meaning that the geometry was first optimized at HF/6-31G* level, and then partial

113

charges were computed from a single-point run, using the CHelpG scheme12 at the B3LYP/6-311++G(2d,2p) ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

114

level of theory; information regarding these molecules is reported in the Supporting Information. This

115

approximation was chosen over the standard OPLS-AA force field calculation (MP2/aug-cc-pVTZ//HF/6-

116

31G*) due to a better stability of DHI+ and DHI+/− when using the 6-311++G(2d,2p) basis set. All the studied

117

models are summarized in Table 1. The number of functional silicate (DHI+ plus DHI+/−) units and structural

118

silicate (SI3 plus SI3−) units was determined from the experimental concentrations of AO-DHI+ iodide and

119

TMOS. It was assumed that the precursors went through complete hydrolysis and fully condensed to DHI+ or

120

SI3 (or its conjugated bases). On the other hand, the ratios of SI3 to SI3− units and DHI+ to DHI+/− units were

121

estimated from a species distribution analysis conducted at pH 9, based on the acidity constant of the silanol

122

group, having in mind that pKa decreases by 1-2 units with high methanol contents.

123

For the systems A and C, the initial state was obtained by inserting into the boxes the respective number of

124

units at random positions using the packmol package13. Regarding the initial state of the B system, we took

125

the final state of the A system and removed a number of water and methanol molecules. This frame, after

126

NpT ensemble equilibration, was set as the initial state of the B system (more details on this process will be

127

given below). Initial box dimensions were estimated considering the molecular weight and the density of

128

each of the components of the mixture. After energy minimization using both steepest-descent and L-BFGS

129

methods included in the GROMACS package, a temperature annealing was performed in the NVT ensemble

130

for 1ns, reaching a temperature of 600 K, so as to ensure a proper mixing and gather three random

131

independent initial configurations. These were, subsequently, used as starting configurations for the three

132

independent MD equilibration runs needed to test the reproducibility of the simulations. Before the

133

production stage, ~100 ns of simulation time in the NpT ensemble were taken to equilibrate the system. As

134

we expected, most of the intermolecular interactions after 100 ns reached a stable value (plots of the energy

135

are reported in the SI, Figure S10), save for the largest groups present in the mixture (e.g.: SI3−, DHI+/− and

136

DHI+). Finally, production runs of 100 ns were performed in the NpT ensemble for data collection.

137

Observable properties were sampled every 2 ps, from which total averages and standard deviations for each

138

run were computed. The equations of motion were integrated using the Verlet leapfrog algorithm14, with a

139

time step of 2 fs. Typically, the temperature (T) was kept fixed at 298 K by applying the velocity rescaling

140

thermostat15, and whenever necessary, the pressure (p) was held constant at 1 bar by using the Parrinello-

141

Rahman scheme16, 17. The time constant used for the Parrinello-Rahman coupling was set to1 ps. Periodic

142

boundary conditions were applied in all three Cartesian directions. For the water molecules, the Transferable

143

Intermolecular Potential four-point model (TIP4P)18 was applied. The non-bonded electrostatic interactions

144

were calculated using a sixth-order Particle Mesh Ewald (PME) method19 beyond a cutoff radius of 0.9 nm.

145

The Lennard−Jones was calculated within a cutoff radius of 0.9 nm with the help of a neighbor list, updated

146

every 10 time steps. A dielectric permittivity, εr, equal to one was used. Statistical and trajectory analysis of

147

the simulations were performed with the programs included in GROMACS, while visualizations were made

148

with VMD20. The analysis consisted essentially in the comparison of the interaction energies between the

149

relevant molecule and ion types, the calculation of radial distribution functions (RDF), coordination numbers

150

(NB), diffusion coefficients (D) and clustering analysis. Both the interaction energy between two groups of


Page 4 of 51

Page 5 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


151

molecules and the energy of interaction within the group were obtained as the simulation average of the sum

152

of short range LJ and Coulomb energies. The RDF between different types of molecules has been calculated

153

as: 〈 〉

= 〈 〉 ,

154

155 156

where refers to the average density of particle B at a distance r, around the particle A, and loc

157

refers to the density of particle B averaged over all spheres around particles A with a maximum radius (rmax)

158

which was half of the box length. The RDF function is additionally averaged on all particles of type A

159

present in the system, and also averaged over the trajectory (simulation time).The g_rdf function included in

160

the GROMACS package calculates the RDF in different ways. The normal method is around a (set of)

161

particle(s), the other methods are around the center of mass of a set of particles or to the closest particle in a

162

set. Here the RDFs were calculated by taking the centers of mass of the particles using the mol_com option

163

of g_rdf function included in the GROMACS package. The coordination numbers (NB) of a particle or atom

164

B around a particle (atom) A were calculated by integrating the radial distribution function gAB(r) between the

165

center of A and the first local minimum, rm:

166

= 4

167 168

where ρB refers to the density of species B (expressed in units of molecules per volume). The cluster analysis

169

was performed using the g_cluster package included in the GROMACS software. This utility can cluster

170

structures using several different methods. We determined structures from the trajectories of the runs using

171

the single linkage which add a structure to a cluster when its distance to any element of the cluster is less

172

than the cutoff. We performed the cluster analysis using cutoff values (i.e., the largest distance to be

173

considered in a cluster) between 0.3 and 1.50 nm. The diffusion coefficient of the mixture components is

174

calculated from the Einstein relation (mean-square displacement, MSD):

175

= 〈| ! − 0| 〉.

176 177 178

where ri is the center of mass positions of the molecules. The MSD is averaged over molecules, and in order

179

to improve the statistics, several restarts r(0) were used along the trajectory.

180 181

3. RESULTS AND DISCUSSIONS

182 183

3.1. Foreword.

184

In this work, we have simulated a rather complex system, comprised of new ionic compoundsSI3−, DHI+/−

185

and DHI+, that had never been modeled before using MD, together with cyclicoxysilane SI3, NAP−, Na+ and



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

186

I− counter ions, as well as the solvents water and methanol. When the complexity and variety of species

187

under consideration is increased, the probability of accumulating a significant amount of calculation errors

188

and uncertainties is also increased. Therefore, it was necessary to assess the validity of the potential

189

considered for the simulations. The main limitation for doing it directly is that the real systems present a

190

huge variety of oxysilane species at a certain moment of the process, making really hard to perform

191

exhaustive atomistic simulations to predict experimental properties of a typical Si based sol-gel process. For

192

this reason, the experimental validation of the model had to be partially indirect. We have chosen OPLS-AA

193

potential including fine-tuned parameters by Price et al.17 to better describe typical sol-gel reagents. First of

194

all, the accuracy of the OPLS-AA potential is extensively verified over small molecules such as water and

195

methanol, so further testing was considered unnecessary for these substances. This potential was previously

196

studied by our group6 for the cyclic SI3 in a mixture which included damascenone, TMOS, water and

197

methanol. In that work, the calculated volumetric mass densities of the different species considered were in

198

good agreement with experimental data, the difference in all the cases being equal to or less than 6.8%

199

(worse case for TMOS).

200

The cationic ORMOSIL derivative, (in the forms DHI+/− and DHI+) is perhaps the most difficult molecule to

201

simulate in our systems due to the presence at the same time of oxysilane rings, dehydroimidazolium cationic

202

moiety and additional permanent charges. The DHI+ is the less investigated structure in our system, and its

203

modeling and validation will be explained more in detail. To compare simulations with experimental data,

204

the ORMOSIL precursor AO-DHI+ iodide was synthesized and characterized. The difference between the

205

synthesized cation AO-DHI+ and DHI+ is the substitution of triethoxysilane and trimethoxysilane by SI3

206

moieties, obtainable after the aforementioned idealized hydrolysis and condensation process. Since both

207

methoxysilane and ethoxysilane parameterizations have been previously validated, a good agreement

208

between experimental and simulated density would provide a satisfactory indication of the adequate

209

parameterization used for the DHI+ substructure. For the simulations, a pure AO-DHI+ iodide system was

210

considered, consisting of 500 ion pairs in a cubic box. The density of the ionic liquid measured in the

211

laboratory was 1180 ± 40 kg /m3, at 298 K, while the density calculated by MD was 1260 ± 1kg /m3. This is

212

a non-negligible difference (6%), but below the 6.8% previously cited. Therefore, we concluded that the

213

model is reasonably good, considering the complexity of the molecule, the potential impurities of the

214

experimental sample (which are not presented in the simulated one) and the inherent limitations of the MD

215

approximation).

216

The calculation of the diffusion is not straightforward because it is very sensitive to modeling conditions

217

such as the force field, time step, and NpT algorithm21. Generally speaking, the analysis of the diffusion and

218

the trajectories of the simulations described in the next sections reveal a clear pre-colloidal character of the

219

mixtures, expressed by the presence of large aggregates of the silicate species (SI3, SI3−, DHI+ and DHI+/-)

220

characterized by a very low mobility (0.004-0.02*10-5 cm2/s) compared to the solvent molecules (0.2-0.9*10-

221

5

222

S7-S9). This behavior was expected and yet reported from other authors22 for the SI3 species alone.

cm2/s). The values of the diffusion coefficients are reported as Supporting Information (Table S5, Figures

223 ACS Paragon Plus Environment

Page 6 of 51

Page 7 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


224

3.2. Model A.

225

As depicted in Figure 1, model A corresponds to the hypothetical intermediate stage reached by the complete

226

hydrolysis of the initial precursors and formation of the SI3 and DHI+ condensates and their alkaline

227

conjugates SI3- and DHI+/-, obeying to the initial experimental concentrations and solvent composition5 (cf.

228

Table 1). As previously referred to, we performed three independent runs for each system. In the case of the

229

A model, the template molecule exhibited an impressive affinity for the DHI+ and a very good one for the

230

DHI+/−. In fact, the analysis of the intermolecular interactions occurring during the simulation clearly shows

231

that major interactions occur between the group pairs DHI+/NAP−, DHI+/SI3− and DHI+/DHI+/−, as reported

232

in Figure 3. As expected, the most unfavorable interaction during the simulation was between the DHI+

233

groups (+2831 kJ/mol), whereas the most favorable interaction energy was recorded for the pairs DHI+/NAP−

234

(−903 kJ/mol), thus indicating a great affinity between this later group. Regarding the other groups, it is

235

important to underline the favorable SI3−/SI3 interaction (−313 kJ/mol) and the unfavorable energy for the

236

pair SI3−/SI3− (+1493 kJ/mol). In addition, is important to underline the good affinity between DHI+ and

237

DHI+/− with SI3− (−826 and −676 kJ/mol respectively). This is a relevant result because it points to an

238

interconnecting trend that leads to the co-polycondensation between the different silicate units in the mixture,

239

thus driving towards the formation of a homogeneous backbone. These results are in agreement with the

240

RDFs. Indeed, analyzing and comparing the graphics of this calculation it is notable the high affinity of the

241

template molecule for the DHI+ and for the DHI+/−, as reported in Figure 4A. The sharp and high peak for the

242

pair NAP−/DHI+ at 0.4 nm clearly shows that the template affinity for the xerogel is optimal. On the other

243

hand, the template molecule shows a very low preference for the counter ion as reported in Figure 3. These

244

results are further confirmed by the NB analysis which clearly shows the preference of the template molecule

245

for the DHI+ and DHI+/− groups, as reported in Table 2. Regarding NB, it is important to underline that this

246

number is calculated using the RDF first global minimum. Thus, it is strictly correlated to the distance in

247

which the minimum is located. In our case, the NAP−/DHI+ pair has a NB mean value of 0.2, but the

248

minimum is at a very close distance (0.5 nm) whereas, for example, for the NAP−/Na+ pair the coordination

249

number is 0.7, but it was calculated at a distance of 1.2 nm, hence highest coordination number does not

250

obviously mean better affinity. Moreover, in all the runs the template does not exhibit appreciable affinity for

251

the SI3 and SI3−, in fact the RDF analysis shows very low peaks compared to the NAP/DHI+ and at a greater

252

distance (0.8nm), as reported in Figure 4D. The SI3 shows a great affinity for itself and for SI3−, in fact the

253

RDF in the Figure 3 reveals a very high and sharp peak for the pair SI3/SI3− at 0.6 nm. In addition, the peaks

254

of the pairs SI3/SI3− and SI3−/SI3− are high remarking the good affinity between these groups. As expected,

255

NB is in agreement with a mean value of ~ 0.6 and ~0.3 for SI3/SI3 and SI3/SI3− pairs, respectively. All these

256

results are in agreement with previous work22 where it was reported the trend of the SI3 to form aggregates;

257

the trend of SI3 and SI3− to form clusters is confirmed by our simulations; the typical cluster of SI3 in this

258

simulation is normally formed by two molecules and along the simulations the number of clusters are usually

259

3 as depicted in Figure 5. Regarding the DHI+ and DHI+/− they show a general trend to form small aggregates

260

of 2 and 4 units (Figure 5 D, E and F), thus indicating that the backbone may potentially grow and the ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

261

building blocks tend to interact, which is expected during the sol-gel process.

262

Overall, the results presented are compatible with a successful imprinting process. In fact, along all the

263

simulation time (100 ns) the template molecules remained near the silica aggregates, mostly those containing

264

the dehydroimidazolium group, which constitutes the most favorable scenario leading to an effective

265

imprinting process. That, of course, does not represent a direct evaluation of the quality of the simulations,

266

but unfortunately the real mixtures are even more complex and it would not be feasible to isolate

267

experimentally any of the calculated parameters. It is only possible at the moment to highlight the

268

consistency between the simulations results and the final imprinting efficiency observed experimentally. In

269

this regard the model B simulations (discussed in the following), both with template and without template

270

inclusion (non-imprinting mixture) leaded us to conditions relatable to nearly gelled mixtures where an

271

additional type of indirect evaluation of the simulations could be made.

272

Another reason for moving to model B (and also C) concerns the very high computational demand (2 ns/day)

273

of model A, since the system was composed of a very large number of units (>8000). In view of even more

274

demanding biphasic future simulations, models B and C would represent significant savings due to removal

275

of ca. 4 fifths of the solvent molecules in model B and reduction of the charged units in model C. In addition,

276

in our opinion, these large systems, in a near future, could be simulated using a coarse graining23

277

methodology and/or with GPU24 based computing, but it is essential to previously simulate the system at an

278

atomistic level in order to gain a better knowledge of the mechanisms regulating these systems.

279 280

3.3. Model B.

281

Experimentally, the mixtures represented by the previous model, may take a very long time to gelify (days or

282

weeks). One common way to speed-up the process is by allowing a slow evaporation of the solvent. The B

283

model was built according to that experimental practice. The same number of units was kept relatively to the

284

A model, except that a number of methanol and water molecules were removed. To roughly estimate the

285

number of water and methanol molecules remaining after this process, the vapor pressure of water (3.2 kPa)

286

and methanol (17.0 kPa) at 298 K, the respective molar fractions in the mixture, and the Raoult law were

287

considered, and also a stepwise batch evaporation model was applied. Basically, a methanol loss of 85% was

288

considered to be achieved in five equalitarian steps. Step 1 consisted on the estimation of how much water

289

would be lost, taking the initial molar fractions and one fifth of the methanol evaporation. In step 2, with the

290

molar fractions, as recalculated after step 1, at two fifth of the methanol evaporation, a second water loss

291

estimation was obtained, and so on.

292

The new model became then much less demanding in terms of computation load (< 1500 units, 20 ns/day),

293

the volume of the simulation boxes being reduced from ~81 to ~47 nm3. Model B could, thus, easily solve

294

the problem of computational effort required by model A. Although the condensed silicate species remain the

295

same, this new model may also be regarded as describing the sol-gel mixture in a later stage of the sol-gel

296

process, which corresponds to a step closer to the gelification point. Here the aggregation level between the

297

silicates would expectedly be much higher, and at the same time the solvent became more polar (more water-

298

rich). Therefore, this model was considered as providing different, but equally useful, information. At a


Page 8 of 51

Page 9 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


299

glance, the removal of solvent has indeed heavily affected the behavior of the system, exhibiting important

300

changes. In general the RDFs confirm the trends reported in the A model, but the three replicas aren’t as

301

similar as the reported for A model. In fact, the template molecule still has better affinity for the DHI+ but

302

shows an increasing affinity for DHI+/− and Na+, as can be seen in Figure 6. Some RDFs differ significantly

303

between each run, but in any case, the template molecule still presents a preferred affinity for the DHI+

304

species confirming the trend shown in the A system. In the case of the first replica (Figure 6A), for the pair

305

NAP−/DHI+, we have a peak at 0.6 nm very close to a maximum at 0.75 nm characterized by a large plateau.

306

In the same figure the RDF of the pair NAP−/Na+ shows a lower peak at 0.6 nm, underlying a growing

307

affinity between these two pairs. This is the most important difference between this system and the A system.

308

For the NAP−/DHI+- pair we recorded the first peak at 0.3 nm, but in this case the affinity was lowered

309

compared with the other groups. In the other two replicas (Figure 6B and 6C) the RDFs of NAP−/DHI+,

310

NAP−/Na+ and NAP−/DHI+/− present a similar behavior. In fact, the RDF of NAP−/DHI+ exhibited two peaks,

311

the first one at 0.4 nm and the second one, very sharp, at 0.6 nm pointing out to a good affinity between these

312

two groups and thus, agreeing with an effective imprinting of the template molecule on the backbone. As in

313

the first replica, the RDF of NAP−/Na+ presents a well-defined peak at 0.65 nm remarking the change in the

314

affinity between the template and its counter ion. Regarding the NAP−/DHI+/−, the first peak with a small

315

plateau at 0.7 nm was recorded, highlighting, as in the other replicas, a low affinity for this pair. These data

316

are in agreement with the coordination number, where the same behavior as in the A system was observed,

317

but with a higher interaction between the template and DHI+/−, as reported in Table 2. Consistently the

318

coordination number for the NAP−/DHI+ and NAP−/DHI+/− pairs is ~ 0.2, denoting a general good affinity

319

between the template and the two species, as previously shown by the RDFs and ~ 0.6 for the NAP−/Na+, as

320

expected considering the RDF analysis. This growing affinity could be explained considering the lower

321

affinity of the template for DHI+. We have to underline here that for these two pairs the coordination number

322

has been calculated at a greater distance than for the NAP−/DHI+, thus it is normal that their values are

323

higher.

324

In agreement with these data, the interaction energies between the template and DHI+/DHI+/− are

325

considerably increased, −1155 and −330 kJ/mol, respectively. In addition, the pairs DHI+/DHI+/−, DHI+/SI3−,

326

DHI+/−/SI3− and DHI+/−/SI3 show a growing interaction energy with respective values of −1428, −925,

327

−1051and −765 kJ/mol, shown in Figure 7. The growing affinity between this species and DHI+/−, SI3 and

328

SI3− means that the backbone continues its trend to grow homogeneously as we should expect to be in a real

329

system. On the other side, the pair DHI+/DHI+ is the group with the most unfavorable energy, +1894 kJ/mol,

330

followed by SI3−/SI3−, +1131kJ/mol, and SI3/SI3 with +964 kJ/mol. Regarding the pair SI3/SI3−, in this case

331

we have recorded the most favorable interaction, −468 kJ/mol, between these molecules which is reflected in

332

the RDF, reported in Figure 8. As it can be seen, there is a low variability between the three replicas. In fact,

333

in the first replica, Figure 8A, a high and sharp peak for each group pair is present in the graph; it is notable

334

that the SI3/DHI+/− pair has a peak at a very close distance, 0.4 nm, whilst for the SI3/SI3 is present at 0.5 nm

335

and for the SI3/SI3− at 0.6 nm, but these later are higher than the former. In the second replica the pair ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

336

groups have the peaks at the same distance but the heights are very different, the SI3/SI3− is the highest

337

whilst, it is notable that the SI3/SI3 is the lowest. In the third replica the SI3/SI3 is the highest whilst, the

338

SI3/DHI+/− is the lowest. These changes are reflected by the coordination number (Table 2) which shows a

339

notable drop of affinity for the pair SI3/SI3 and a slight drop for the pair SI3/SI3−. Consistently, the cluster

340

analysis shows a higher propensity of SI3 and SI3− to form large aggregates confirming the general trend of

341

these units as reported in Figure 9A, B, and C. The cluster analysis also reveals that DHI+ and DHI+/− in this

342

system are forming very large aggregates (Figure 9D, E and F).This fact is really relevant because it endorses

343

the other result we obtained highlighting that this simulation is capable of reproducing the behavior of this

344

sol-gel system at a stage closer to gelification where a large and homogenous backbone is presented. This,

345

together with the affinity of the template for the DHI monomers, is really indicative of an effective

346

imprinting effect of the NAP− on the backbone. All these results are supported by the coordination number of

347

the pair SI3/SI3 with a mean value of ~0.4 for the SI3/SI3 and 0.2 for the SI3/SI3−. Regarding the DHI

348

monomers, the NB confirms the precedent results with a notable increase compared with the A model; in fact,

349

we have recorded a value of 0.4 for the pair DHI+/DHI+ and 0.5 for the pair DHI+/DHI+/−, as can be seen in

350

Table 2.

351

Given the relatively low computational demand of model B, we decided to simulate also a variant of this

352

model (model B-NI) consisting on the simple removal of the template (NAP-), representing thus a model of

353

the mixture for the preparation of the non-imprinted xerogel, a reference material routinely prepared

354

experimentally for ascertaining the imprinting effect. The comparative analysis of the interaction energies

355

(Table S6) and the side by side presentation of snapshots taken from the trajectories of the simulations of

356

models B and B-NI (Figure S6) were quite elucidative of the effect of the presence of the template. The

357

snapshots clearly evidence the structuring effect of the template upon the silicate aggregates, substantiated by

358

a dispersive action. The presence of the template, bearing good affinity with the gel backbone, implied a real

359

interposition between the network-forming units, leading to a much less compact structure able to

360

accommodate the template. The values of the interaction energies corroborated such dispersive modulating

361

effect by the template. Indeed, globally it was observed that the presence of the template impeded a higher

362

level of interaction between the network-forming units, as that occurring within the non-templated mixture.

363

Such textural effect of the naproxen template must translate into the porosity features of the final xerogel,

364

and, accordingly, the BET N2 adsorption analysis of the experimental counterparts of models B and B-NI

365

showed5 that, although both materials presented very small pore volume (imprinted 0.0038 mL/g vs. non-

366

imprinted 0.0013mL/g) the imprinted material presented significantly higher porosity.

367

The results obtained with these models are thus compatible with the behavior of mixtures that led to

368

successfully imprinted xerogels, indicating that this kind of simulation can be useful in the study of such

369

systems, at least in qualitative terms. The template vs. non-templated mixtures comparison showed the

370

expected templating effect, in agreement with the RDF analysis and NB computation. Thus, in our opinion

371

this kind of approach has the potential to be used to prepare models of nearly gelified imprinted gels that

372

could be used to evaluate the recognitive properties of a MIP polymer against a NIP.

373 ACS Paragon Plus Environment

Page 10 of 51

Page 11 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


374

3.4. Model C.

375

The last simulated model was built to compare the behavior of the system in absence/reduction of some ionic

376

species in the system, lowering the computational demand required by A model, in this case by reducing the

377

number of charged species. This system was built using the same number of molecules of water and

378

methanol as in the A model, but did not comprise the SI3- and the DHI+/− species. In other words, we wanted

379

to ascertain if the acid-base equilibrium of silanol could be suppressed without implicating erroneous

380

conclusions or lose some important information. In order to respect the same stoichiometry of the A model,

381

the system comprised 20 molecules of DHI+ and 18 of SI3, as reported in Table 1. As for the other systems,

382

we performed the analysis of the intermolecular interactions. Again, the pair NAP−/DHI+ had the most

383

favorable interaction energy, −1337 kJ/mol, clearly showing the good interaction for this pair. In addition, it

384

is remarkable the good affinity between the DHI+/SI3 pair with a favorable interaction energy of −627

385

kJ/mol. On the other hand, the pair DHI+/DHI+ is the group with the most unfavorable interaction energy

386

followed by the pair SI3/SI3, +3638 and +1584 kJ/mol respectively, as can be seen in Figure 10. Another

387

important consideration when comparing the interaction energies of the three simulated models is that, in this

388

case, the standard deviation along the three simulations is the lowest, due to the reduced number of charged

389

species present in the box. Consistently with the energies, the RDF analysis reveals that the NAP− still has a

390

good affinity for the DHI+, but at the same time has a better association for its counter ion present in the

391

system (Na+) as reported in the Figure 11A. This change could not be totally unexpected because in this

392

system the DHI+ cannot easily interact with its anionic counterpart (DHI+/−) due to the absence of the anionic

393

charge. Thus, in this system the silicate species form a smaller and less homogeneous backbone that is more

394

disperse. Hence, the Na+ and NAP− could better diffuse thorough the solvent and interact between each other.

395

In the case of DHI+ the simulation shows a drop affinity for the template and a growing affinity for the SI3,

396

as reported in Figure 11B. These changes are induced by the growth affinity of the NAP− for the Na+.

397

Furthermore, analyzing the trajectories and the cluster analysis the trend of the SI3 to form large aggregates

398

is confirmed; in fact the cluster analysis along the three replicas shows a clear trend of the SI3 to form

399

clusters formed by 4 or 5 monomers as reported in Figure 12. This result is endorsed by the RDF analysis

400

which reveals that the SI3 has a great affinity for itself and an increased affinity for the DHI+ (Figure 11B).

401

Taking all these data, it is remarkable that the template molecule still has a good affinity for the DHI+,

402

highlighting also in this simulation the favorable interactions towards an effective imprinting of the template

403

molecule on the silica backbone. In qualitative terms, this model allowed to reach similar trends for the most

404

relevant interactions (template-template, template-functional silicates, template-structural silicates and

405

functional silicates – structural silicates) concerned with the imprinting process. However, the gain in

406

computational performance in absence of DHI+/− and SI3− was not that significant (2 ns/day for system A vs.4

407

ns/day for system B). In addition and more important, the results of this model compared to the A system

408

clearly reveal that an explicit representation of the ionic species is essential to study and understand the

409

behavior of such complex system. In fact, due to the absence of the DHI+/− in this system, it was impossible

410

to evaluate if a polycondensation trend occurring between major species (DHI+, DHI+/−, SI3 and SI3−) is

411

present. In this context, it is preferable to simulate a system with a major number of species, in order to gain ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

412

also an understanding on the mechanisms that may involve all the possibly charged species. With these

413

results in mind, we thought it was not necessary to simulate the B system with the same approximation we

414

did from system A to C, since there would not be a significant speed-up of the simulations and, on the other

415

hand, an important loss of information may occur.

416 417

4. CONCLUSIONS

418

This paper represents a significant step forward for future studies of the molecular imprinting process in a

419

sol-gel system, through atomistic MD simulations. It is important to underline that MD simulations can

420

provide molecular level details and uncover the atomistic basis of the interaction between the xerogel and the

421

template molecule. The simulated systems were rather complex and included a new functional silicate in its

422

cationic and zwitterionic form, DHI+ and DHI+/−, and the anionic form, SI3−, of the well-studied SI3. In this

423

context, the main objective of this paper was focused on the study of the affinity of DHI+/DHI+/− for the S-

424

Naproxen, the template molecule. For the first time, these new silicates have been simulated in a series of

425

atomistic MD simulations with an explicit representation of all the ionic species accordingly to the pH of the

426

system; in this sense, the physical plausibility of the parameterization of new molecules (to the OPLS-AA

427

force field) was successfully checked in order to assess the reliability of the systems and simulations. The

428

analysis of the simulations by means of RDF, NB, and the energy of interaction of the groups, have

429

demonstrated the high affinity of the xerogel to the template molecule even after loss of solvent (model B).

430

In this sense, the A models were the most representative in providing the clear imprinting effect of the

431

template molecule on the xerogel; thus, all the results obtained clearly show that the S-Naproxen has a great

432

affinity for the DHI+ and a good affinity for the DHI+/−. In the B models, with a much higher density of

433

silicate units and other conditions resembling a nearly gelified sol, the imprinting effect was not so marked,

434

but the association between the NAP− and DHI+ was always the highest compared with other groups. This

435

model allowed us also to observe a clear texturing dispersive effect of the template on the silicate network in

436

other to be accommodated within that network. That was no less than testifying the own process of molecular

437

imprinting taking place. Finally, it is important to notice that, although the results of the simulations could

438

not be validated directly against any property measured experimentally, the scenarios observed by MD,

439

especially for model B, were totally consistent with the behavior observed for the materials prepared in the

440

lab, namely a good imprinting effect ascertained by an imprinted/non-imprinted adsorption and porosity

441

comparison.

442 443

ASSOCIATED CONTENT

444

Supporting Information

445

Atomic point charges and geometry used for the DHI+, DHI+/− and SI3− molecules. This material is available

446

free of charge via the Internet at http://pubs.acs.org

447 448

ACKNOWLEDGMENTS


Page 12 of 51

Page 13 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


449

This work received financial support from the European Union (FEDER funds through COMPETE) and

450

National Funds (FCT, Fundação para a Ciência e a Tecnologia) through projects Pest-C/EQB/LA0006/2013

451

(REQUIMTE) and Pest-C/QUI/UI0081/2013 (CIQ). The work also received financial support from the Eu-

452

ropean Union (FEDER funds) under the framework of QREN through Project NORTE-07-0124-FEDER-

453

000067-NANOCHEMISTRY. RC acknowledges also FCT and the European Social Fund for financial sup-

454

port (Grant SFRH/BPD/80605/2011). To all financing sources the authors are greatly indebted.

455 456 457 458 459 460 461 462 463 464

References

465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494

(1). Collinson, M. M. Sol-Gel Strategies for the Preparation of Selective Materials for Chemical Analysis. Crit. Rev. Anal. Chem. 1999, 29 (4), 289-311. (2). Zhang, H.; Ye, L.; Mosbach, K. Non-covalent molecular imprinting with emphasis on its application in separation and drug development. J. Mol. Recognit. 2006, 19 (4), 248-59. (3). Díaz-García, M. E.; Laínño, R. B. Molecular Imprinting in Sol-Gel Materials: Recent Developments and Applications. Microchim Acta 2005, 149 (1-2), 19-36. (4). Vasapollo, G.; Sole, R. D.; Mergola, L.; Lazzoi, M. R.; Scardino, A.; Scorrano, S.; Mele, G. Molecularly Imprinted Polymers: Present and Future Prospective. International Journal of Molecular Sciences 2011, 12 (9), 5908-5945. (5). Kadhirvel, P.; Azenha, M.; Shinde, S.; Schillinger, E.; Gomes, P.; Sellergren, B.; Silva, A. F. Imidazolium-based functional monomers for the imprinting of the anti-inflammatory drug naproxen: Comparison of acrylic and sol–gel approaches. J. Chromatogr. A 2013, 1314 (0), 115123. (6). Azenha, M.; Szefczyk, B.; Loureiro, D.; Kathirvel, P.; Cordeiro, M. N.; Fernando-Silva, A. Molecular dynamics simulations of pregelification mixtures for the production of imprinted xerogels. Langmuir 2011, 27 (8), 5062-70. (7). Azenha, M.; Szefczyk, B.; Loureiro, D.; Kathirvel, P.; Cordeiro, M. N.; Fernando-Silva, A. Computational and experimental study of the effect of PEG in the preparation of damascenoneimprinted xerogels. Langmuir 2013, 29 (6), 2024-32. (8). Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26 (16), 1701-1718. (9). Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics 1983, 79 (2), 926-935. (10). Azenha, M.; Szefczyk, B.; Loureiro, D.; Kathirvel, P.; Cordeiro, M. N. l. D. S.; FernandoSilva, A. n. Molecular Dynamics Simulations of Pregelification Mixtures for the Production of Imprinted Xerogels. Langmuir 2011, 27 (8), 5062-5070. (11). M. J. Frisch, G. W. T., H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537

Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc., Wallingford CT, . Gaussian 09. 2009. (12). Breneman, C. M.; Wiberg, K. B. Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis. J. Comput. Chem. 1990, 11 (3), 361-373. (13). Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M. PACKMOL: a package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30 (13), 215764. (14). Hockney, R. W.; Goel, S. P.; Eastwood, J. W. Quiet high-resolution computer models of a plasma. J. Comput. Phys. 1974, 14 (2), 148-158. (15). Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. The Journal of Chemical Physics 2007, 126 (1), -. (16). Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52 (12), 7182-7190. (17). Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50 (5), 1055-1076. (18). Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225-11236. (19). Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. The Journal of Chemical Physics 1995, 103 (19), 8577-8593. (20). Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 14 (1), 33-38. (21). Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D. Molecular Dynamics Simulation of Liquid H2O, MeOH, EtOH, Si(OMe)4, and Si(OEt)4, as a Function of Temperature and Pressure. The Journal of Physical Chemistry A 2001, 105 (10), 1909-1925. (22). Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D. Molecular Dynamics Simulation of Methanolic and Ethanolic Silica-Based Sol−Gel Solutions at Ambient Temperature and Pressure. The Journal of Physical Chemistry A 2001, 106 (1), 130-148. (23). Gohlke, H.; Thorpe, M. F. A natural coarse graining for simulating large biomolecular motion. Biophys. J. 2006, 91 (6), 2115-20. (24). Ruymgaart, A. P.; Elber, R. Revisiting Molecular Dynamics on a CPU/GPU system: Water Kernel and SHAKE Parallelization. Journal of chemical theory and computation 2012, 8 (11), 4624-4636.

538 539 540 541 542 543 544 545 546 547 ACS Paragon Plus Environment

Page 14 of 51

Page 15 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571

572 573 574 575 576 577 578


Tables

Table 1. Description of the Simulated Systems

model

methanol:water ratio

NAP-

A B C

25:1 5:1 25:1

10 10 10

COMPONENTS (NUMBER OF UNITS) DHI+ DHI+/- SI3 SI3- wate methanol r 10 10 20

10 10 -

9 9 18

9 9 -

310 230 310

Table 2. Coordination numbers calculated for selected groups in the simulation. A B -

NAP /NAP NAP-/DHI+ NAP-/DHI+/NAP-/SI3 NAP-/SI3NAP-/Na+ SI3/DHI+ SI3/ DHI+/SI3/SI3 SI3/SI3DHI+/DHI+ DHI+/DHI+/-

0.07 0.20 0.11 0.04 0.01 0.74 0.11 0.11 0.59 0.33 0.16 0.10

0.14 0.21 0.26 0.66 0.05 0.65 0.09 0.16 0.43 0.24 0.38 0.49

579 580 ACS Paragon Plus Environment

7736 1130 7736

Na+

I-

29 29 10

20 20 20

C 0.20 0.10 0.09 0.05 0.11 0.79 -


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600

601 602 603 604 605 606 607 608 609 610 611 612

Figures

Figure 1

Figure 1. Scheme of the formation of naproxen-imprinted xerogel obtained from a mixture of TMOS and AO-DHI+ iodide. The intermediate state shows the main species selected for the simulations, obtained after an idealized process of complete hydrolysis and cyclic trimerization. In the case of Model A the conjugated bases of the two silicate trimerswere also included (not shown) and the amount of solvent (methanol:water 25:1) reflected the initial experimental conditions.In model B we kept the number and nature of thesilicatedspecies and template as in model A, but the amount of solvent was much reduced (to ca. 1 fifth) and its proportion changed to methanol:water5:1). It was intended here to study mixtures with characteristics of silicate density and solvent composition typical of a stage closer to the gelification. Model C was equivalent to model A with the difference that we did not include the conjugated bases of the shown cyclic trimers.

Figure 2


Page 16 of 51

Page 17 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

613 614 615 616 617

618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638


Figure 2. Molecular structure of DHI+ (A), DHI+/- (B), SI3 (C), SI3- (D), NAP- (F).

Figure 3

Figure 3. Interaction energies for the most important groups obtained in the SI3DHISI3:25:1 model.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658

Figure 4


Page 18 of 51

Page 19 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

659 660 661 662 663 664 665 666 667 668 669 670 671 672 673


Figure 4. RDF calculated from the trajectories of the three replicas ofSI3DHISI3:25:1 run. A, B and C are representative of the three replicas for group pairs NAP-/DHI+, NAP-/Na+, NAP-/DHI+/-.C, D and E are representative of the three replicas for group pairs NAP-/SI3, NAP-/SI3-.

Figure 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694

Figure 5. Cluster analysis of the SI3, NAP and SI3- calculated from the trajectories of the SI3DHISI3:25:1 runs, each image is representative of one replica.

Figure 6


Page 20 of 51

Page 21 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

695 696 697


Figure 6. RDF calculated from the trajectories of the SI3DHISI3:25:1 run for the atom-type pairs, NAP/DHI+, NAP/-DHI+, NAP/NA.

Figure 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739


Figure 8


Page 22 of 51

Page 23 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

740 741 742


Figure 8. RDF calculated from the trajectories of the SI3DHISI3:25:1 run for the atom-type pairs SI3/SI3, SI3/SI-, SI3/DHI+.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759

Figure 9. Cluster analysis of the SI3 calculated from the trajectories of the SI3DHISI3:5:1 runs, each image is representative of one replica.

Figure 10


Page 24 of 51

Page 25 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801


Figure 10. Interaction energies for the most important groups obtained in the SI3:25:1 model.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

802 803 804 805

Figure11. RDF calculated from the SI3:25:1 run.


Page 26 of 51

Page 27 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834


Figure 12. Cluster analysis of the SI3 molecule; 1, 2 and 3 are the numbers of the replica.

TABLE OF CONTENTS (TOC) ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

835

836 837


Page 28 of 51

Page 29 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Scheme of the formation of naproxen-imprinted xerogel obtained from a mixture of TMOS and AODHI+ iodide. The intermediate state shows the main species selected for the simulations, obtained after an idealized process of complete hydrolysis and cyclic trimerization. In the case of Model A the conjugated bases of the two silicate trimerswere also included (not shown) and the amount of solvent (methanol:water 25:1) reflected the initial experimental conditions.In model B we kept the number and nature of thesilicatedspecies and template as in model A, but the amount of solvent was much reduced (to ca. 1 fifth) and its proportion changed to methanol:water5:1). It was intended here to study mixtures with characteristics of silicate density and solvent composition typical of a stage closer to the gelification. Model C was equivalent to model A with the difference that we did not include the conjugated bases of the shown cyclic trimers. 369x242mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Molecular structure of DHI+ (A), DHI+/- (B), SI3 (C), SI3- (D), NAP- (F). 796x322mm (96 x 96 DPI)


Page 30 of 51

Page 31 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Interaction energies for the most important groups obtained in the SI3DHISI3:25:1 model. 581x250mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. RDF calculated from the trajectories of the three replicas ofSI3DHISI3:25:1 run. A, B and C are representative of the three replicas for group pairs NAP-/DHI+, NAP-/Na+, NAP-/DHI+/-.C, D and E are representative of the three replicas for group pairs NAP-/SI3, NAP-/SI3-. 755x874mm (96 x 96 DPI)


Page 32 of 51

Page 33 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Cluster analysis of the SI3, NAP and SI3- calculated from the trajectories of the SI3DHISI3:25:1 runs, each image is representative of one replica.

826x800mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6. RDF calculated from the trajectories of the SI3DHISI3:25:1 run for the atom-type pairs, NAP/DHI+, NAP/-DHI+, NAP/NA. 391x824mm (96 x 96 DPI)


Page 34 of 51

Page 35 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



576x238mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8. RDF calculated from the trajectories of the SI3DHISI3:25:1 run for the atom-type pairs SI3/SI3, SI3/SI-, SI3/DHI+.

420x815mm (96 x 96 DPI)


Page 36 of 51

Page 37 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9. Cluster analysis of the SI3 calculated from the trajectories of the SI3DHISI3:5:1 runs, each image is representative of one replica.

758x819mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 10. Interaction energies for the most important groups obtained in the SI3:25:1 model. 702x244mm (96 x 96 DPI)


Page 38 of 51

Page 39 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure11. RDF calculated from the SI3:25:1 run. 334x495mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 12. Cluster analysis of the SI3 molecule; 1, 2 and 3 are the numbers of the replica. 283x216mm (96 x 96 DPI)


Page 40 of 51

Page 41 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure S1. DHI+ molecule 309x184mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure S2. DHI+/- molecule 335x168mm (96 x 96 DPI)


Page 42 of 51

Page 43 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure S3. NAP molecule 261x155mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure S4. SI3 molecule 206x191mm (96 x 96 DPI)


Page 44 of 51

Page 45 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure S5. SI3- molecule 231x201mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure S6. Snapshot of the three replicas of the models B and B-NI 1163x843mm (96 x 96 DPI)


Page 46 of 51

Page 47 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure S7. Graph of the diffusion in the Model A 1816x1204mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure S8. Graph of the diffusion in the Model B 1900x1246mm (96 x 96 DPI)


Page 48 of 51

Page 49 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure S9. Graph of the diffusion in the Model C 1930x1333mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure S10. Graphs of the Kinetic energy, Lennard-Jones short range, Coulomb short range, Potential energy and Total energy of the three replicas of the Model B. 1483x905mm (96 x 96 DPI)


Page 50 of 51

Page 51 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


346x259mm (96 x 96 DPI)


Molecular Dynamics Simulations of Complex Mixtures Aimed at the

Recommend Documents