Automated optimization of water-water interaction parameters for a coarse-grained model.

We have developed an automated parameter optimization software framework (ParOpt) that implements the Nelder-Mead simplex algorithm and applied it to a coarse-grained polarizable water model. The model employs a tabulated, modified Morse potential with decoupled short- and long-range interactions incorporating four water molecules per interaction site. Polarizability is introduced by the addition of a harmonic angle term defined among three charged points within each bead. The target function for parameter optimization was based on the experimental density, surface tension, electric field permittivity, and diffusion coefficient. The model was validated by comparison of statistical quantities with experimental observation. We found very good performance of the optimization procedure and good agreement of the model with experiment.

Introduction

Water plays a fundamental role in a wide range of fields of study, including geology, atmospheric sciences, chemistry, biology, and physics.[1] Because of its relevance and ubiquity, water is the most common solvent in experimental and computational studies of biological systems.[2] The hydrophobic effect, the segregation of molecules based on relative interactions with water, governs membrane self-assembly, protein folding, and many other biological processes.[3] Despite its seemingly simple structure, the water molecule forms intricate and dynamic networks of hydrogen bonds that give it unique bulk and interfacial properties.[4,5] The wealth of anomalous characteristics that water exhibits relative to most other liquids underscore its uniqueness. In part because of its enigmatic character, it has been the subject of many theoretical and modeling efforts in the past few decades.[6] Reproducing these anomalous characteristics poses a serious challenge for any effort to attain a detailed understanding of water dynamics and structure.

Quantum mechanical models of water based on first principles (ab initio) precipitate the evolution of the electronic structure along with nuclear coordinates over time. The resolution of these methods is at the level of molecular orbitals. These ab initio models produce radial distribution functions (RDFs) and rotational/vibrational spectra that strongly agree with experimental data.[7,8] The computational complexity of quantum models restricts systems to a range too limited to determine most statistical bulk quantities. These properties, such as the diffusion coefficient and thermodynamic variables, are more accurately modeled by semiempirical approaches. Water models with atomic-level detail are the most common for simulation of biomolecules and their aggregates. To reproduce a wide range of complex water behavior, several models that consist of four (Bernal–Fowler,[9] TIP4P,[10] and TIPS2[11]) and five (TIP5P[12]) interaction sites have been developed. In the interest of reducing computational cost, most simulations use models (SPC,[13] SPC/E,[14] and TIP3P[15]) that are composed of three partially charged atoms connected with rigid bonds and a Lennard-Jones (6-12) intermolecular potential between oxygen atoms. Each model is reasonably accurate within its domain of applicability, but none reproduce a large portion of the anomalous properties of water. Various models[16] include flexible bonds that allow the water dipole moment to change with the surrounding environment.

Coarse-grained (CG) models are able to achieve longer time scales by averaging over rapid inter- and intramolecular motions, thus permitting a larger time step. CG models are able to effectively represent a larger number of molecules by reducing the number of interaction sites considered. Because of the unusual properties of water, the development of a CG model that reliably describes the bulk and medium-range properties is a challenging project. The development of models of molecular assemblies with lower complexity that reproduce important qualities of water[17,18] is ongoing. For an extensive review, see Coarse-Graining of Condensed Phase and Biomolecular Systems.[19] Johnson et al.[20] have identified limitations of CG models by demonstrating the lack of transferability of a CG pair potential across different states. The work also shows that the pair potential cannot simultaneously resolve all the properties of the reference system for a given state. Despite these limitations, there is extensive development of various CG models for water.[21−32] One common CG protocol is to group a number of water molecules into a single bead with its center as the interaction site.[21−24,26,27,32] This class of CG models carries no charge. Electrostatic interactions are implicitly incorporated into the effective pairwise interaction potential which is either a Lennard-Jones (LJ)[21−24] or Morse-like function.[27] The free parameters in pair potentials are trained against thermodynamic properties such as density, surface tension, and solvation free energy. Molinero and Moore[32] built a coarse-grained water model around the tetrahedral properties shared by water, silicon, and carbon by adapting the Stillinger–Webber potential originally developed for silicon. The model maps one water molecule onto one interacting bead and describes important properties of water across a wide range of temperatures. Using a clustering algorithm, Hadley and McCabe derive the effective potential through fitting to structural data.[17] The models of Shelley et al.[21] and Chiu et al.[27] have large isothermal compressibility. Using a LJ 12-4 potential, Shinoda et al.[24] have designed a CG water that has density and surface tension comparable to those of water. Electrostatic interactions can be introduced into CG water models with multiple charged sites per bead.[25,28−31] These models can describe the dielectric properties of water. Despite the incorporation of electrostatic interactions, both the BMW and polarizable MARTINI models predict density different from experimental measurements.[28,29] The CG representation of Darré et al.[30] associates approximately 11 water molecules with four tetrahedrally interconnected beads (WT4). The model explicitly accounts for long-range electrostatics. The parameter set of the WT4 model is tuned to match the experimental density. On the other hand, the isothermal compressibility and surface tension are modeled less accurately.[30] The recent GROMOS CG[31] is a 5-to-1 mapped water model with two electrostatic interaction sites. It delivers good results for density, surface tension, and dielectric constant as compared to those of real water but yields a coefficient of thermal expansion 1 order of magnitude higher than the experimental value. Adaptive multiscale models allow transfer between simulation domains with differing granularity, from the quantum level up to CG.[33−35] Adaptive methods face unique challenges because of the variability in the degrees of freedom over the span of a simulation. In continuation of our previous CG development based on the Morse potential,[27] we upgrade our previous CG water model (termed CSJ) with polarizability and a Morse-like potential with more flexibility in the landscape of pairwise interaction parameters. We will refer to the new model as modified Morse coarse-grain (MMCG).

Model

A CG mapping of four water molecules to one CG water bead was used for the current model. For interactions between water beads, we have further refined our previous model,[27] which used the standard Morse potential, to meet two goals: (1) to improve compressibility without sacrificing other features of the potential and (2) to build into the interaction potential a smoothing function that permits fairly accurate simulations using a shorter cutoff. The choice of pairwise potential gives flexibility in the functional form without introducing any extra simulation costs. The coarse-grained potential must model the average of several different interatomic and intermolecular forces. Therefore, the functional form is solely based on a phenomenological understanding of the system rather than a first-principle derivation, such as the r–6 behavior of the dispersion force included in the Lennard-Jones potential. A modified Morse potential (eq 1) was used to describe interactions between water bead centers.wherewhere R and De are the potential well location and depth, respectively. Compared to the Morse potential, the MMCG potential replaces α, the scaling parameter in the exponent, with two separate distance-dependent variables, α(r) and β(r). This form decouples the repulsive (eq 1a) and attractive (eq 1b) parts of the potential. The shape of the repulsive term as r tends to 0 is determined by α0. The exponential factor tends toward αL as r approaches R. Parameters βR and βc provide similar behavior for the attractive part of the potential. Beyond the cutoff distance rc, the potential is taken to be 0. At this writing, we have not conducted an exhaustive search for the value for rc that optimally combines accuracy and efficiency, but the conventional value of 1.6 nm gives good results with reasonable efficiency. Comparing systems representing equal numbers of water molecules, the MMCG model achieves a >1 order of magnitude improvement in performance over an atomistic system. The polarizability of CG water beads is modeled like that of the polarizable MARTINI water model,[28] through the addition of a harmonic angle potential term, with spring constant K0 and equilibrium angle θ0, among three charged points. This model of polarizability allows for not only dipole but also higher moments of the electrostatic energy. Figure 1 provides an illustration of the MMCG model. The central point, which acts as the interaction site for the modified Morse potential, carries negative charge Q, while the two outer points each carry charge −Q/2. A mass equal to that of four water molecules is distributed evenly among the three points. The distance between the outer points and the central point is a tunable parameter I. The change in the angle (θ) between the points represents the polarizable nature of the bead. The set of potential parameters {De,R,α0,αL,βR,βc,θ0,K0,Q,I} determines the space on which optimization will be performed.

Figure 1

Illustration of the MMCG model. Four water molecules map onto one CG bead. Each CG bead has three charge sites. The central site with charge Q provides the interaction site for the modified Morse potential. Outer charge sites (with charge −Q/2) are connected to the central site by a fixed bond of length I. The central angle is a determined by a harmonic potential with equilibrium angle θ0.

Methods and Results

Parameter Optimization

Parameter optimization is exploration of a parameter space to determine extreme values of a target function. In the case of molecular dynamics force field optimization, the set of potential parameters defines the space and comparison of model predictions with either ab initio data or experimental measurements defines the target function to be minimized. We have developed a new software package, Parameter Optimizer (ParOpt), for general optimization problems (ParOpt is available for download under the GNU public license at https://csmlabfs1.cas.usf.edu/Sites). The software provides various optimization algorithms, including an orthogonal steepest descent and the Nelder–Mead simplex-based method. See Figure 2 for a schematic representation of ParOpt. Target functions defined for molecular dynamics force field optimization present a challenging problem. Sources of difficulty include high dimensionality, discontinuities, possible stochasticity due to random numerical error and finite sample size effects, and difficulty in defining a derivative. Thus, in the context of molecular dynamics parameter optimization methods such as simplex-based, approaches with the least demanding conditions on the target function are the most suitable. Nelder–Mead, a commonly used simplex-based optimization algorithm, iteratively evolves M + 1 points on the M-dimensional space using four basic moves.[36]

Figure 2

Schematic representation of the ParOpt optimization framework. ParOpt provides an environment in which different optimization methods can be accessed by defining the generation of the target function.

Reflect

Replace the highest point with a point reflected about the centroid of the remaining points.

Expand

If the reflected point is lower in value than all points in the centroid, consider a point further from the centroid.

Contract

If the reflected point would remain the highest, consider a point nearer to the simplex.

Shrink

Move all points nearer to the lowest point.The M + 1 individual points in the M-dimensional space are denoted by P, while subscripts h and l denote the lowest and highest points, respectively. The centroid is defined byValues used for optimization of water parameters are those suggested by Nelder–Mead:[23] α = 1.0, γ = 2.0, β = 0.5, and δ = 0.5. Figure 3 graphically illustrates the Nelder–Mead simplex transformations.

Figure 3

Schematic illustration of Nelder–Mead steps. Arrows denote points that will be replaced and the new point to be considered: (a) the reflection of the highest point about the centroid of the remaining points, (b) the further expansion of the reflection move, (c) the contraction along the line joining the centroid and the highest point, and (d) the shrinking of the simplex by shifting all points nearer to the lowest point.

The Nelder–Mead algorithm replaces the highest point in the simplex with one with a lower target function value by applying these moves. A Nelder–Mead step begins with calculation of a reflect move. If the target function value for this new reflected point is lower than all points in the simplex, an expansion move is attempted. If this expanded point is lower than all points in the simplex, it is accepted. If the reflected point is not higher than the highest point in the new simplex, it is accepted and replaces the highest point in the simplex. If the newly accepted reflected point is still the highest point in the simplex, a contraction step is attempted. If the contracted point is higher than the highest point in the simplex, then a shrink move is applied, replacing all points except the lowest. If the contracted point is lower than at least one point, it is accepted into the simplex. This process is repeated until the stopping criterion is met. The stopping criterion for the Nelder–Mead method is a cutoff on the root-mean-square target function value over the entire simplex. The number of occurrences of each move in an optimization depends on the characteristics of both the simplex and the target function and is therefore problem specific. The parameter space for the MMCG water model is a 10-dimensional bounded space. Minimal and maximal values (see Table 1) for each parameter were chosen heuristically to restrict target function evaluations to reasonable combinations of parameters. The original description of the Nelder–Mead algorithm has no set method for boundary conditions on the parameter space, though suggestions are made for modification of functional forms or assignment of large “penalty” numbers for points outside of the boundary. We chose to replace points outside the boundary with the nearest point within the allowable space.

Table 1

Optimal Parameters for a Set of Final Parameters from Nelder–Mead Optimization Compared with Results from the Previous Model (CSJ)

parameter	CSJ value[27]	optimized value (this work)	minimal value	maximal value
De (kJ mol–1)	3.4	3.742	2.0	5.0
R (nm)	0.629	0.560	0.55	0.61
α0	7	15.545	7.0	20.0
αL	7	13.323	7.0	20.0
βR	7	10.976	7.0	20.0
βc	7	3.198	1.0	7.0
θ0	–	131.406	105.0	145.0
K0 (kJ mol–1)	–	93.918	20.0	180.0
Q (e)	–	–1.126	–1.6	–1.06
I (nm)	–	0.141	0.1	0.25

For optimization of the MMCG water model, the target function was defined by the weighted percent error in comparison of CG simulation results with four experimental quantities: density, dielectric constant, diffusion coefficient, and surface tension. To ensure that all comparison data contribute to the target function and thus the evolution of the simplex, percent errors are assigned scalar weights to yield similar orders of magnitude among all quantities. With this goal in mind, density was weighted 100 times higher than surface tension, while permittivity and diffusion coefficient were weighted 10 times higher. Experimental values used in training are listed in Table 2.

Table 2

Training Data Showing Results for the Optimized Point Compared with Experimental Target Data

property	experiment	training result
ρ (g cm–3)	0.996[39]	0.993
D (×10–9 m s–2)	2.597[40]	2.61
ε	76.8[41]	76.24
γ (mN m–1)	71.2[42]	78.43

Density

Density was computed by taking the ratio of the mass of the system and the ensemble average volume.

Diffusion Coefficient

The standard practice in comparing diffusion coefficients between CG models and experiment is to assume that diffusion for N-sized clusters scales as 1/N; i.e., D = D1/N, where D is the diffusion coefficient for an N-sized cluster. This assumption neglects any interatomic correlation. The GROMOS CG work[31] compares the CG bead diffusion coefficient with the diffusion coefficient for the center of mass of clusters in the SPC water model. In that work, static clusters in the atomistic representation were produced by adding distance restraints between oxygen atoms, which produces diffusion coefficient scaling close to the inverse behavior usually assumed. In the work presented here, clusters are constructed similarly, though unrestrained. To calculate diffusion coefficient scaling, clusters were constructed from each molecule and its nearest N – 1 neighbors. The MSD for such clusters was calculated using a window averaging method with eq 7. Though the method presented here may not be generalizable to all types of clusters, it provides an improvement over current assumptions for homogeneous systems. Figure 4 shows the diffusion coefficient scaling for SPC/E water. The plot for atomistic water shows significant deviation from the scaling behavior of uncorrelated clusters, especially at larger values of N. Therefore, approximations leading to 1/N scaling are more appropriate at smaller cluster sizes. At a bead size of four water molecules, the difference is significant. On the basis of this analysis, we use a scale factor (s) of 3.16 to compare the CG diffusion coefficient with experiment, instead of the usual factor of 4.where x(t) is the position of the central bead at time t and the broken brackets denote an ensemble average.

Figure 4

MSD scaling for SPC/E water. The line shows the usual practice of taking the diffusion of an N cluster beginning N times slower than the single molecule. SPC/E was determined by simulation of SPC/E at 303 K. Clusters were chosen by the N – 1 nearest neighbors for each water molecule, and the center of mass diffusion was calculated.

Relative Permittivity

The relative permittivity (ε) was calculated using a Clausius–Mossotti-like equation with a reaction field with an infinite dielectric constant (conducting boundary conditions).[37]where M is the total system dipole, ε0 is the permittivity of free space, ⟨V⟩ is the ensemble average system volume, T is the target temperature for the thermostat, and ε0 is the permittivity of the vacuum.

Surface Tension

The surface tension was taken from the GROMACS internal calculation:where the z-axis is normal to the interface, P terms are the pressure tensor diagonal elements, Ns is the number of interfaces, and L is the size of the system along the z-axis.

The optimization procedure involves many target function evaluations. These evaluations, at a point in parameter space, involve the execution of two MD simulations. This places practical limits on the size and length of simulations used for optimization. Therefore, systems used for determination of the diffusion coefficient, density, and electric field permittivity contained 512 CG beads, corresponding to 2048 water molecules. This system is denoted as S1. Because of large fluctuations, the calculation of surface tension requires a larger system: 4096 CG beads, implying 16384 water molecules (system S2). System S2 was constructed, at every evaluation of the target function, with a slab near the average density computed from S1 in contact with a vacuum.[38] System S2 was simulated under the constant number, volume, and temperature (NVT) ensemble. More details on simulation methods are found in Appendix A.

The first step in optimization is the construction of an initial simplex. In this work, optimization started with a simplex consisting of random points within the domain of parameter space that satisfied constraints. This simplex was optimized until convergence. To further improve the accuracy of the final point, the optimization was restarted with an initial simplex consisting of the converged point and random points around it. This procedure was then iterated two more times, with the final optimization resulting in the point presented in Table 1.

Figure 5 shows the mean and root-mean-square (rms) target function value versus the Nelder–Mead step for the final optimization iteration. Because of the roughness of the hypersurface generated by the target function, the simplex may occasionally include a vertex with an abnormally high value, as observed in the mean and rms of the second shrink step of Figure 5. Figure 6 shows the evolution of parameter values over the final run, while Figure 7 shows the change in computed physical properties. The figures demonstrate the convergence of both the data and parameter values. Further, Figure 6 demonstrates that although parameters did encounter the boundaries over the course of the optimization, the final simplex converged away from the boundary. Thus, the parameter value boundaries chosen were not unreasonable.

Figure 5

Optimization results. The mean error, rms error, and rms cutoff as a function of Nelder–Mead step are shown. Reflection (□), contraction (○), and shrink (■) steps are labeled at the corresponding mean error point.

Figure 6

Optimization results. Parameter trajectory. Variation in parameters as a function of Nelder–Mead step. Dashed lines indicate minimal and maximal values.

Figure 7

Optimization results. Training data [density (ρ), electric field permittivity (ε), diffusion coefficient (D), and surface tension (γ)] (—) as a function of Nelder–Mead step compared with target values (---).

The final set of parameters was chosen as the point in the final converged simplex with the lowest target function value (see Table 1). The decoupling of repulsive and attractive parameters, the added distance dependence of the exponential factors, and the added electrostatic interaction have allowed for a decrease in the equilibrium distance, while maintaining the experimental density. The large differences between exponential parameters in the new and original models hint that the extensions have added important flexibility. The simulation results corresponding to the final set of parameters are listed in Table 2.

Validation

When validating the results of an optimization procedure, one cannot consider the training data alone. The suitability of the force field will depend on matching data that the model was not explicitly trained to reproduce. Therefore, validation simulations were performed to verify the reliability of the optimized parameters. Because validation is not performed iteratively, it is not subject to the same size constraints as optimization. All systems for validation consisted of 110592 CG beads, which corresponds to 442368 water molecules. See Appendix A for more details on simulation methods.

The density, diffusion coefficient, relative permittivity, and surface tension were determined from the larger validation systems using the same methods that were used during optimization [see section 3.1; in addition, bulk thermodynamic quantities were also calculated (see Table 3)].

Table 3

Validation Data Comparing Experimental Data, an Existing CG Model (polarizable MARTINI), an Atomistic Model (SPC/E), and the Previous Model (CSJ) to This Work (MMCG)a

property	experiment	polarizable MARTINI[28]	SPC/E	CSJ[27]	MMCG (this work)
ρ (g cm–3)	0.996[39]	1.043	0.998[14]	0.998	0.993 ± 0.0005
D (×10–9 m s–2)	2.597[40]	2.5	2.5[14]	4.3	3.07 ± 0.03
ε	76.8[41]	75.6	70.7 (at 298 K)[49]		74.17 ± 0.06
γ (mN m–1)	71.2[42]	30.5		71.0	78.73 ± 0.2
α (×104 K–1)	3.21[50]				6.89 ± 1.3
κT (×10–6 bar–1) (NPT)	44.75[51]		34.07[44]	170.0	68.95 ± 10
κT (×10–6 bar–1) (NVT)	44.75[51]		41.41[44]		57.16 ± 0.5
Hvap (kJ mol–1)	44.0[52]			38.4	48.52

Quantities computed are density ρ, diffusion coefficient D, dielectric permittivity ε, surface tension γ, coefficient of thermal expansion α, isothermal compressibility κT, and enthalpy of vaporization Hvap. Values computed both during optimization and for validation differ because of the size differences between simulated systems.

Coefficient of Thermal Expansion

The coefficient of thermal expansion, α, was calculated using the finite-difference method.[43]Two constant number, pressure, and temperature (NPT) ensemble simulations at temperatures of 303.15 K (T1) and 308.0 K (T2) were performed. The ensemble averages of the resulting densities were used for ρ1 and ρ2, respectively.

Isothermal Compressibility

Isothermal compressibility κT was calculated using two methods:[44] (i) a finite difference method for NVT simulations at different densities (ρ)and (ii) volume fluctuations in an NPT ensemble simulationwhere P is the pressure on system i and ΔV2 is the variance of the volume. For eq 11, the following densities were used: ρ1 = 0.993 g cm–3, and ρ2 = 0.963 g cm–3.

Enthalpy of Vaporization

The enthalpy of vaporization was calculated from the interaction energy between CG beads (VInter) and the intermolecular interactions that occur within the bead VS of 128.5 kJ mol–1 determined from quantum mechanical calculations of the water tetramer binding energy.[27,45]where the factor of 4 is due to the level of CG.

The average dipole moment of an individual bead was calculated for the MMCG water model. Because a CG bead represents four water molecules, the total dipole moment for an atomistically detailed SPC/E water molecule along with its three nearest neighbors was computed for comparison. The MMCG model yields an average dipole moment of 4.1 ± 0.25 D, whereas the four-water cluster dipole moment from the SPC/E water model was 5.3 ± 1.0 D; density function theory (DFT) calculations for an isolated cluster produce no dipole moment.[46] The effective polarizability (α0) for a CG bead was calculated using eq 14.[47] The MMCG model yields a value of 63 au compared with a value of 39.33 au given by DFT for an isolated four-water cluster.[48] The disparity in these values is most likely a result of differences between clusters in the bulk and in isolation.

Conclusion

To automate force field optimization, we have developed a general optimization framework ParOpt. The software utilizes various methods to locate minima of the defined target function on parameter space. Parameters are systematically varied to determine the most accurate set of parameters. We have applied this general framework to optimize the parameters of a polarizable CG water model with nonbonded interactions given by a modified Morse potential via the Nelder–Mead algorithm. The model has more flexibility in the functional form than most CG models, allowing it to match a wide range of experimental properties (see Table 3).

Any model in a classical molecular dynamics framework is a simplification of a physical system and therefore involves a reduction in the degrees of freedom in the description of the interactions. Any such modeling attempt would be nonunique. Though the set of parameters that govern the model are mathematically independent, their effects on observable data and therefore the target function are not necessarily uncorrelated. Therefore, because of this interdependence of the parameters, the optimization is not necessarily underdetermined despite an apparent abundance of free parameters to match the target data. The question of uniqueness is nontrivial to answer and is in fact a matter of ongoing research (e.g., in the context of economic models[53]). For complex systems in general, there is no straightforward, computationally tractable, and general solution to this problem; rather, heuristic approaches must be developed for the particular system being studied. In light of these factors, it may be the case that our parameter set presented in this paper is one of many equally good sets. The point we have chosen was a suitable point that matched the training data well. Additionally, it can be seen from the good agreement with the validation data, which is independent of the training data, that the resulting force field parameters are well-tuned.

There are fundamental limitations to the explanatory power of CG water models. RDFs are notoriously difficult for to match. An artifact common to most CG models is the long-range correlation between CG water beads not seen in experiment or atomistic models. Figure 8 shows that the current model suffers from that same limitation. Likely as a result of the high degree of interbead correlation at long distances, the model has the propensity to spontaneously solidify at room temperature, a weakness found in other CG water models.[23] This phenomenon would not be expected in a heterogeneous system, because long-range correlations would be broken by the presence of other interactions. Additionally, spontaneous freezing can be avoided by periodically assigning random velocities chosen from a Maxwell distribution, a method similar to one suggested by Harvey et al. to prevent the overpopulation of low-frequency modes of a molecular dynamics system.[54] Because dynamic and structural properties of water clusters vary with temperature, the internal states of CG beads would, as well. Thus, CG models cannot necessarily be assumed to transfer to different temperatures. We believe that the interaction parameters for CG models would be more accurately treated as being temperature-dependent. The set of potential parameters has been tuned for one temperature and is not expected to accurately simulate behavior far from the target temperature or in thermodynamic ensembles without constant temperature control.

Figure 8

RDF. Comparison between g(r) for the atomistic SPC/E water model (—) and the MMCG water model (---). The CG model has greater long-range coordination than the atomistic model.

The predictions of the model agree well with experiment for structural, dynamic, and bulk properties. The diffusion constant shows that the model has accurate dynamics. The dielectric constant indicates the validity of the parameters that determine electrostatics and polarizability. Density ensures the correct spatial scale of the CG beads. Accurate surface tension indicates an accurate representation of the strength of the interaction between beads. The thermodynamic quantities predicted by the model are in good agreement with those of experiment, validating the fitness of the parameters.