extfieldsofstudy:
Medicine
Computer Science
Physics
Chemistry

text:
Development of effective stochastic potential method using random matrix theory for efficient conformational sampling of semiconductor nanoparticles at non-zero temperatures

In this work, the development and implementation of the effective stochastic potential (ESP) method is presented to perform efficient conformational sampling of molecules. The overarching goal of this work is to alleviate the computational bottleneck associated with performing a large number of electronic structure calculations required for conformational sampling. We introduce the concept of a deformation potential and demonstrate its existence by the proof-by-construction approach. A statistical description of the fluctuations in the deformation potential due to non-zero temperature was obtained using infinite-order moment expansion of the distribution. The formal mathematical definition of the ESP was derived using functional minimization approach to match the infinite-order moment expansion for the deformation potential. Practical implementation of the ESP was obtained using the random-matrix theory method. The developed method was applied to two proof-of-concept calculations of the distribution of HOMO-LUMO gap in the water molecule and solvated CdSe clusters at 300K. The need for large sample size to obtain statistically meaningful results was demonstrated by performing $10^5 $ ESP calculations. The results from these prototype calculations demonstrated the efficacy of the ESP method for performing efficient conformational sampling. We envision that the fundamental nature of this work will not only extend our knowledge of chemical systems at non-zero temperatures but will also generate new insights for innovative technological applications.

I. INTRODUCTION
The goal of this work is to develop an efficient theoretical and computational method for performing ensembleaveraged properties of semiconductor nanoparticles (NPs). Semiconductor NPs such as quantum dots and rods have important technological applications. Because of the quantum confinement effect [1], their optoelectronic properties can be tuned by adjusting parameters such as the size, shape, material composition, and ligands of the NP. Consequently, NPs have promising applications in photovoltaics, [2][3][4][5][6] light emitting devices, [7][8][9] charge and energy transfer processes [10][11][12][13], and biological labeling. [14,15] One of the important environmental factors that influence their optical and electronic properties is the temperature, which has been the focus of numerous theoretical and experimental studies. Early experimental investigation performed by Joshi and co-workers in 2006 looked at this temperature dependence of band gap energies for colloidal CdSe/ZnS core/shell nanocrystals of different diameters. They found that over a temperature range of 10K to 300K, the band gap energies decreased with increasing temperature for all nanocrystal sizes [16]. Recent experimental results on the temperature-effect on quantum dot (QD) optical properties have also been obtained by Savchenko and co-workers in 2017. They studied the first exciton peak shift as a function of temperature in InP/ZnS core/shell QDs in the range of 6.5K to 296K. Similarly, to Joshi et al. they found that the peak both broadens and shifts toward larger wavelengths with increasing temperature [17]. In 2011, Chen and co-workers performed nonadiabatic molecular dynamics simulations on both spherical and elongated CdSe quantum dots using time-domain Kohn-Sham theory. Their results confirm the inverse temperature dependence on band gap energy in quantum dots. They also found that the temperature dependence of the relaxation rate of excited electrons was higher in elongated quantum dots [18]. In the same year, Fischer et al. performed time-dependent density functional theory (TD-DFT) calculations on a Si 29 H 24 quantum dot trajectory. They found that the band gap energy of the QD shifted significantly throughout the dynamics simulation [19]. More recent quantum dynamics simulations utilizing TD-DFT and the tight-binding approximation performed by High and co-workers looked at the excited state energy transfer in a porphyrin lightharvesting device. Their simulation showed that it is essential to couple the nuclear and electronic degrees of freedom to accurately simulate excited state energy transfer [20]. Liu and Jakubikova have also investigated the electron transfer rate in a pyridine-TiO 2 light harvesting assembly, and found that the transfer rate was faster when sampling from thermally accessible structures at room temperature compared to a single 0K optimized structure [21].
In colloidal systems, the NPs are surrounded by solvent molecules which provide a constant-temperature heatbath for the NP. Constant-temperature ensembles such as NVT or NPT are used extensively for calculation and analysis of equilibrium properties. In molecular dynamics simulations, various methods are used to implement the thermostat, such as Berendsen [22], Langevin [23], and Nosé-Hoover [24,25]. Investigation of temperature dependence is complicated and computationally expensive because it requires electronic structure calculations on a large number of deformed structures. This imposes serious limitations on the size of the systems that are being investigated.
Disorder arising from temperature or solvent effects can also impact the observed properties of organic molecules and systems as well. Much theoretical and computational work has been performed in this area. In 2001, Kwasniewski and coworkers calculated UV-vis absorption spectra on various frames throughout molecular dynamics simulations of trans-stilbene using the ZINDO-CIS method. They found that thermal motion had a strong impact on the HOMO to LUMO transition energy. At temperatures less than 400K, the energy levels displayed pronounced Gaussian thermal broadening of the related band, up to 30 nm [26]. In 2003, Lewis and co-workers performed MD simulations combined with tight-binding density functional theory to calculate the electronic states of a model DNA double helix as the molecule underwent room temperature thermal fluctuations. They found that the population density of the HOMO state over time changed dramatically due to these thermal fluctuations [27]. In 2006, Morth and Autschbach investigated the temperature dependence of the optical rotation of fenchone using TD-DFT. They found that optical rotation in this system increases linearly with temperature over the range of 273K-373K, and concluded that the temperature-dependent vibrational effects were largely responsible for the temperature-dependent rotational effects in rigid organic molecules [28].
Many researchers have also investigated these solvent and temperature effects on organic chemical systems experimentally. In 2002, Ariu et al. studied the temperature dependence on photoluminescence quantum yields (PLQYs) of poly(9, 9-dioctylfluorene) films. They found that the PLQYs increased as the temperature of many of these systems were reduced. In 2005, Dauphas and co-workers investigated β-casein protein aggregation and properties as a function of concentration, temperature, and calcium content, characterized by dynamic light scattering, confocal microscopy, and fluorescence spectroscopy. Over a temperature range of 283K to 323K, they saw that the fluorescence absorption profiles of the protein were very different when the concentration and temperature were changed. In 2014, Mani et al. looked at the fluorescent properties of the salt of 2,6-diaminopyridinium with dihydrogen phosphate. They found that as a result of proton transfer from phosphoric acid in solution to the pyridine nitrogen, the band gap between the HOMO and LUMO states is smaller than compared to 2,6-diaminopyridine. This finding was confirmed by a red shift that was observed in its absorption spectrum [29].
In this work, we present the effective stochastic potential (ESP) method to address the computational bottleneck associated with performing electronic structure calculations on a large number of molecular structures needed for obtaining ensemble-averaged quantities. The central idea of the ESP method is to approximate the explicit interactions between the system and the bath degrees of freedom by a statistically equivalent effective potential. In the ESP method, the construction of the effective potential is performed using the random-matrix theory (RMT) method. Before we describe this new method, a brief background in the history and usage of RMT will be provided for context. Originally, RMT was designed by Wigner in 1951 to deal with the statistics of eigenvalues and eigenfunctions of complex many-body quantum nuclear systems [30]. Wigner later used RMT to investigate the statistical fluctuations of scattering processes. [31] After Wigner, RMT was successfully applied to the description of spectral fluctuations of atomic nuclei [32], complex atoms [33], and complex molecules [34]. The nature of the general theory greatly lends itself to the investigation of highly complex and chaotic systems. Other, more recent examples of how RMT has been applied to physical systems include equilibrium and transport properties of disordered and classically chaotic quantum systems [35], transitions between classical and quantum distributions of hydrogen energy levels [36], and spectral resonances in quartz [37]. Much work has also been done using RMT to investigate the quantum mechanical phenomenon of so-called "persistent current," the perpetual flow of interacting electrons in a static magnetic field [38][39][40][41][42].
There are many forms of RMT that have been used, but in general, RMT is the process of replacing the Hamiltonian of a system with an ensemble of Hamiltonians containing random matrix elements. Most popularly used in RMT is the Gaussian orthogonal ensemble (GOE), where the Hamiltonian is partitioned into two parts: a reference component and a stochastic component where the matrix elements are uncorrelated, normally distributed random numbers. These two components are then coupled to form the total Hamiltonian. Due to their highly complex nature, QDs are a natural candidate for investigation using RMT. One of the earliest investigations of QDs using RMT was done by Jalabert, Stone, and Alhassid in 1992. They developed a statistical theory of the effect of Coulomb blockade in QDs represented as a particle in a finite potential well using random matrices to determine the Coulomb blockade amplitudes. They were able to obtain good agreement with experiment. [43] In 1996,Seba et al. used a random matrix model to describe the conductance of a model QD connected to an ideal lead. To accurately match experimental results of QDs attached to two leads, they introduced a coupling constant to their RMT Hamiltonian. [44] In that same year, Sivan and co-workers performed an experimental investigation into the disorder of ground state energies of populations of GaAs QDs. [45]. More recently, in 2008 Kaplan et al. investigated the fluctuation of two-body screened Coulomb interaction in ballistic QDs using the random wave model (a model related to RMT). They were able to derive analytical expressions for these fluctuations in terms of the linear size of the QD. Despite the Gaussian random single electron basis used in their theory, they found that both two-body and one-body matrix elements followed non-Gaussian distributions. [46] Shankar in 2006, and later in 2008, has discussed a technique for studying QDs using the renormalization group method in conjunction with RMT. [47,48] In 2017, Akimov utilized stochastic Hamiltonians to investigate excited state relaxation in nonadiabatic condensed matter systems [49].
The remaining part of this article is organized as follows. The main theoretical development of the ESP method is described in section II. The derivation includes the definition of the quantum mechanical Hamiltonian (section II A) under investigation, the definition of the deformation potential (section II B), and the conceptual construction of the effective stochastic potential (section II C). The computational details for implementation of the ESP method are presented in section III. Specifically, implementation details on the choice of sampling schemes(section III A), choice of one-particle basis functions (section III B), and stochastic sampling using Gaussian random matrices (section III C) are presented . The developed ESP method was used to perform proof-of-concept calculations of ensemble-averaged HOMO-LUMO gap energies in chemical systems, and results for H 2 O and Cd 20 Se 19 are presented in section IV and the discussion of the results are presented in section V.

A. System setup and definitions
We start by introducing key concepts and definitions that will be used repeatedly in this work. We will start by defining the effective one-particle Hamiltonian as The one-particle Hamiltonian is the starting point for treating electron-electron correlation in many-electron systems.
In the above expression, v eff is the effective one-particle operator for treating electron-electron interaction, which can be approximated using HF, DFT, pseudopotential, MBPT, or some model potential For wave function based methods, the eigenfunctions of the one-particle Hamiltonian provide the one-particle basis for performing, MP2, CI, and, CCSD calculations.
In case of DFT, these one-particle basis functions are used to construct the one-particle density. The derivation presented here is very general and does not make any a priori assumption about the form of the effective potential. Typically, v eff is not known in advance and has to be constructed iteratively using a self-consistent procedure. However, because of the SCF procedure, the final SCF-converged v eff has a functional dependence on the nuclear geometry R Using the SCF converged effective potential for any nuclear geometry R η , we define the following one-particle Hamiltonian We also define a reference structure R 0 and the effective Hamiltonian corresponding to the reference structure is given by f 0 : Typically, the reference structure is selected to be the minimum energy structure. However, the derivation does not impose this as a requirement.

B. Definition of deformation potential
We define a structure R η to be deformed if is different from the reference structure R 0 , and cannot be generated from the reference structure by simple rotational and translational transformation on R 0 . The difference between R η and R 0 is a consequence of deformation in the internal coordinates of the molecular system. Associated with each deformed structure, we define the deformation potential, which is defined as Using the above definition, the deformation potential can be interpreted as the potential which must be added to the reference Hamiltonian to generate the deformed Hamiltonian. The existence of the deformation potential is an exact condition, and can be proved using the following argument. Because f η and f 0 exist, therefore, using Eq. (7) v η def must also exist. Up until this point, we have used abstract Dirac representation for representing the reference and deformed Hamiltonians. However, now we seek to represent the deformed potential as a matrix in a single-particle basis. The matrix representation of the above equation is given as where the matrix elements are defined using as set of single-particle basis functions{φ The choice of the single-particle basis functions is an important one, and in principle, may or may not depend on the nuclear geometry R η . The choice of the single-particle basis in which the deformed potential is represented is an important topic and will be discussed in section III.

C. Conceptual construction of Effective Stochastic Potential
Calculation of the deformation potential becomes computationally expensive for a large set of structures. The effective stochastic potential method is designed to address this computational bottleneck by replacing the exact deformation potential by the effective potential that shares common identical statistical metrics with the deformation potential for a set of structures. Conceptually, the ESP potential can be defined by the following three steps. First, we define a set of all possible unique structures that a chemical system can exist in by {R η } ∞ 0 . Second, for each of the structures, the deformation potential is determined, where the set of all deformed potentials is defined as For each deformed structure R η , we associate a corresponding probability of p η for existence of that structure in a physical system. As the definition of probability implies, p η should satisfy the following two mathematical relationships: and η∈S def p η = 1.
For the remainder of the derivation we will use the compact notation of v η ij to refer to the individual matrix elements of the deformation potential For each matrix element of V η def the mean µ or the expectation value of the distribution is defined as Here, we have used both s (1) ij and µ ij to refer to the mean for convenience. Analogously, the variance σ of the distribution is defined as Using the mean µ ij , we define the m th central moment of the distribution by s m which is evaluated using the following expression: Evaluation of all the moments for each matrix element of the set of deformation potentials implies the following relationship: Using these statistical metrics, we are now in the position to define the effective stochastic potential in matrix representation. We start by defining a symmetric stochastic matrix whose matrix elements are random numbers Z drawn from a probability-distribution function f pdf ij (z): The probability distribution function must satisfy the following two conditions: The probability of finding Z in the interval [a, b] is given by We define the effective stochastic potential as the stochastic potential whose central-moments are closest to the centralmoments of the deformed potential. Mathematically, the ESP potential is defined using the following constrained functional minimization: where the constraints on f pdf ij are given in Equation 20.

III. COMPUTATIONAL AND IMPLEMENTATION DETAILS
In this section, we present the implementation details of the ESP method by introducing additional approximations needed for practical implementation of the method. Practical implementation of the ESP method requires us to work with finite sets of data. As a consequence, only a finite set of structures {R η } M 0 were used for evaluation of the deformation potential and the set of deformation potentials was also finite. Computational implementation also required us to limit the set of moments (s (1) , . . . , s (N ) ) on the set of deformation potentials to a finite number. The various methods of sampling, choice of single-particle basis functions, and determination of the probability distribution function for ESP are presented below.

A. Canonical ensemble sampling using Monte Carlo and molecular dynamics
The probability p η associated with the structure R η depends on the thermodynamic conditions of the physical system, and can be chosen to describe both equilibrium and non-equilibrium conditions. In this work, we are interested in systems at thermal equilibrium, which are well-described using the constraints of canonical ensemble with constant composition, temperature, and volume (N,V,T). For canonical ensemble, the probability is the well-known Boltzmann expression as shown below: where, β = (k B T ) −1 and E η is the energy of the deformed structure and E 0 is the energy associated with the minimum-energy structure. The first step in the determination of the deformation potential is the generation of the set of structures R η . In principle, the set of structures can be achieved using either Monte Carlo (MC) by molecular dynamics (MD) simulations. The key point in both of these approaches is to generate a statistically meaningful sample of thermally accessible structures. In the MC procedure, this is achieved by randomly distorting the optimized structure, calculating its energy, and and calculating the mean and central moments using the Boltzmann-weighted procedure described in Equation 23. In the MD procedure, equilibrium molecular dynamics calculations are performed at fixed N,V,T and the structures are then randomly selected from the equilibrium distribution. It is important to note that Boltzmann weights are only needed for MC sampling, and should not be included in the case of MD sampling. The relevant equations for MD sampling are presented in Equation 24 In the above equation, the superscript (MD) implies that the set of structures were obtained from the MD simulation. B. One-particle basis for representing the deformation potential An important feature of the ESP method is the choice of the one-particle basis functions that are used for representing the deformation potential. In this work, we will use the eigenvectors of the reference Hamiltonian for representing both the deformed and the ESP potential. We start by defining the AO-basis representation of the reference Hamiltonian matrix In the first step, we perform canonical orthogonalization In the next step, we obtain the unitary matrix that diagonalizes the reference transformed Fock matrixF 0 Using U 0 , we transformed all the deformed Fock matricesF η into the eigenbasis of the reference Fock matrixF 0 It is important to note that we are transforming using U 0 instead of U η . Finally, we obtain the matrix representation of the deformed potential in the eigenbasis of the reference Fock matrix using the following expression which is equivalent to subtracting the eigenvalues of the reference Fock matrix from the deformed Fock matrix.

C. Stochastic sampling using Gaussian random matrices
In this work, we have restricted the sampling to only the mean and the variance The above sampling condition is computationally implemented using symmetric Gaussian random matrices. Each element of the ESP matrix is is a Gaussian random number Z as shown in Equation 33 The function f Gauss ij is the Gaussian or normal probability distribution function with mean µ and variance σ, and has the following form The corresponding Fock matrix associated with ESP is given bỹ IV. RESULTS

A. Statistical distribution of HOMO-LUMO gap of H2O at 300 K
We begin with a proof-of-concept calculation for the ESP method, with the goal of demonstrating that the properties of a distribution of thermally-accessible geometries of a chemical system can be accurately reproduced with a stochastic potential derived from that system. We use a single isolated water molecule as our test case. This serves as a useful example of a simple, well-studied polyatomic system. Starting with NVT ensemble, our goal is to calculate the frequency distribution of the HOMO-LUMO gap of water at 300K. The calculations were performed using the canonical ensemble with constant N, V, and T at 300K. The sampling of the structures for this ensemble was performed using the Monte Carlo procedure. A set of 1000 structures of water molecules were generated by introducing random distortions to the bond-distances and bond-angles of the optimized structure of water. For each of the 1000 structures, Hartree-Fock (HF) calculations were performed using cc-pVDZ basis functions, using the GAMESS software package. [50] The probability associated with each of the 1000 structures were calculated using the Boltzmann expression as shown below where E 0 is the energy of the minimum-energy structure. The eigenvalues of the Fock matrix from the converged HF calculations were used to calculate the HOMO-LUMO gap for the 1000 structures.
Using the thermal probabilities p i , the thermally-weighted frequency distribution of the HOMO-LUMO gaps for 1000 structures obtained from the HF/cc-pVDZ calculation at 300K is presented in Figure 1. The mean and variance of the distribution is presented in Table I.

Construction of the effective stochastic potential
The ESP for water at 300K was constructed from a different set of 1000 Fock matrices obtained from the electronic structure calculations. Using Equation 25, all the Fock matrices were first transformed into the eigenbasis of the Fock matrix associated with the minimum-energy structure F 0 . Then, using Equation 31, the set of 1000 deformation potential matrices were obtained. The mean and the variance associated with each matrix element of the deformation potential matrices were calculated using Equation 14 and Equation 15. The means and variances were used in the next step to calculated the matrix elements of the ESP using Gaussian random numbers as show in Equation 34. The ESP Fock matrix was constructed from the ESP matrix using Equation 35. A set of 1000 ESP Fock matrices were generated by stochastic sampling, and the eigenvalues of the ESP Fock matrices were used to calculate the ESP HOMO-LUMO gap

Comparison of the ESP and electronic structure results
The comparison of the mean, the variance and higher central moments of the frequency distributions of HOMO-LUMO gap energies obtained using both Hartree Fock and ESP methods are presented in Table I. Using the mean µ, the p th root of the central moment is defined as It was found that both the mean HOMO-LUMO gap obtained from ESP were in good agreement with the HF results. The results demonstrate that the ESP can successfully capture the fluctuations present in the deformation potential at non-zero temperatures. However, comparison of higher-order moments (Table I) reveal that the two distributions have similar but not identical characteristics. In order to compare the statistics of the Hartree Fock and ESP method distributions, we used Z-score (also called standard score) [51]. Z-score is defined mathematically as where x is a measured value, µ is the mean of the value over a distribution, and σ 2 is the corresponding standard deviation. Z-score is a unitless quantity that is relative to the distribution's standard deviation. A positive Z-score reflects a value greater than the mean, and negative Z-score reflects a value smaller than the mean. A Z-score of 1 corresponds to a value that is exactly one standard deviation greater than the mean of the distribution. Although, the ESP and HF distribution were found to differ in certain regions (Figure 1), the overall features of the both the distributions were found to in good agreement with each other.

Effect of sample size on frequency distribution
The effect of sample size on the ESP results were investigated by performing 10 3 , 10 4 , and 10 5 uncorrelated sets of ESP calculations, and the results are summarized in Table I and Figure 2 . In all cases, the mean and higher-order moments exhibited convergence with the size of the sample. The frequency distributions in Figure 2 show that the noise in the distribution for N sample = 1000 is attenuated by increasing the sampling size to 10 4 .

B. HOMO-LUMO gap of Cd20Se19 in water at 300K
We now turn the discussion to the application of the ESP method to a small CdSe QD (1.2 nm in diameter) in aqueous media at 300K. These systems are excellent candidates on which to apply the ESP method, due to the potentially useful applications of QDs, the complexity of their electronic states, and the impact of solvent environment and non-zero temperature of its optical and electronic properties. and ESP method using increasing number of sampling points.

System setup and construction of the effective stochastic potential
Equilibrium molecular dyanamics simulations with constant N,V,T were performed for the Cd 20 Se 19 cluster with explicit water molecules at 300K. The force-field developed by Rabani et al. [52] was used for the CdSe cluster, and the explicit water molecules were treated using TIP3 force. The MD calculations and analysis of the simulations were performed by the GROMACS molecular dynamics program. The initial structure for the MD calculations was obtained from the minimum-energy structure obtained from electronic structure theory at HF/LANL2Z ECP level. The initial structure was then re-optimized using the CdSe force-field in the MD program in vacuum. The simulation cell was constructed by solvating the CdSe cluster in a water bath consisting of 11763 water molecules. The entire system was equilibrated at 300K by performing MD simulation for 500 ns. The equilibration conditions were verified by monitoring the total energy, total potential energy, and the temperature of the system.
After equilibration, a set of 50 structures was randomly selected from the MD trajectories for construction of the ESP matrix. As a first step, HF calculation using LANL2Z ECP basis functions were performed for the set of structures, and the deformation potential matrices were calculated. The mean and variance of each matrix element were calculated, and the ESP matrix and the ESP Fock matrix were constructed using the procedure described in section III.

Predictions from the ESP calculations
The average HOMO-LUMO gap (µ) and the higher-order moments (Equation 39) obtained using the ESP are presented in Table II. The results from Table II highlight the need for large sampling size to obtain statically meaningful quantities. Specifically, we find that a sample size of N sample = 50 is inadequate for describing the distribution of the HOMO-LUMO gap for both HF and ESP calculations. The results from ESP show that sample size in excess of 1000 are needed to converge the mean and higher-order moments. Comparison between the HF and ESP results show that mean HOMO-LUMO gap is significantly over-estimated in HF set of calculations because of inadequate sample size. The distribution of the HOMO-LUMO gap from the ESP calculations are presented in Figure 3. The distribution was found to be highly asymmetric with respect to the origin with a significant positive bias towards the -Z direction and a long tail in the +Z direction. The frequency distribution from the ESP data ( Figure 3) implies that at 300K, the HOMO-LUMO gap for majority of the structures are lower than the mean gap. The long tail implies, that a very small fraction of the population exhibits very large HOMO-LUMO gaps. The distribution at 300K is compared with the the HOMO-LUMO gap of the minimum-energy structure (vertical black line in Figure 3). The comparison show that the maxima of the frequency distribution is significantly different from the minimum-energy structure. This implies  The vertical bar at Z=2.2 represents the HOMO-LUMO gap energy of the 0K optimized structure with ∆ǫ = 2.7eV. [53] .
that the minimum-energy structure is a poor representative descriptor for the set of thermally accessible structures at 300K. From Figure 3 we conclude that the heating the system from 0K to 300K should decrease the HOMO-LUMO gap of CdSe clusters. This observation from the ESP calculation also is supported by previous experimental [16,17] and theoretical [18,19] studies on CdSe and other semiconductor clusters and quantum dots.

V. DISCUSSION
Sampling constant temperature ensembles: This work demonstrates that the ESP method can be used with either Monte Carlo or molecular dynamics simulations. Specifically, the water calculations were performed using Monte Carlo sampling and the CdSe cluster calculations were performed using classical molecular dynamics simulation. In our experience, classical MD simulation provides a an efficient computational route for sampling NVT and NPT ensembles for solvated molecular systems which explicit solvent molecules. We believe that the study of CdSe cluster in water are particularly relevant to the study of semiconductor nanoparticles, because it has generated new insight into the electronic structure of clusters at non-zero temperatures in aqueous media. The results from the ESP calculations also highlighted the requirement of a large sample set for obtaining statistically meaningful results.
Future directions: Although the proof-of-concept calculations presented here were restricted to HF calculations, the ESP method is can be combined other methods for treating electron-electron correlation. The ESP method generates a manifold of single-particle states which can be used as a starting point for combining ESP with existing many-body theories such as DFT, TDDFT, and Green's function methods. In addition, the ESP method can also be applied to non-equilibrium or non-Boltzmann distributions by appropriate choice of the set of probabilities {p η } (either discrete or continuous) that govern the population of molecular structures for a given chemical environment.

VI. CONCLUSIONS
This work presented the development of the effective stochastic potential method (ESP) as a route to perform efficient conformational sampling for molecular systems at non-zero temperatures. The central idea of the ESP method is to replace the fluctuations experienced by the molecular system by a stochastic potential that provides an equivalent statistical description of those fluctuations. The concept of deformation potential was introduced and the existence theorem for such a potential was presented. The formal mathematical constructions of the ESP potential was presented using a functional minimization approach to capture the fluctuations in the deformation potential. The computational implementation of the ESP was achieved using random matrix theory. Two proof-of-concept calculations were performed using the developed ESP method, by calculating the distributions of HOMO-LUMO gaps of water and solvated CdSe cluster at 300 K. The results from these calculations demonstrate the efficacy of the ESP method in providing an efficient conformational sampling of molecular systems at non-zero temperatures.


extfieldsofstudy: Medicine Chemistry

s2fieldsofstudy: Chemistry Biology

text:
Retention Mechanism of Proteins in Hydroxyapatite Chromatography – Multimodal Interaction Based Protein Separations: A Model Study

Background: Retention mechanism of proteins in hydroxyapatite chromatography (HAC) was investigated by linear gradient elution experiments (LGE). Materials and Methods: Several mobile phase (buffer) solution strategies and solutes were evaluated in order to probe the relative contributions of two adsorption sites of hydroxyapatite (HA) particles, C-site due to Ca (metal affinity) and P-site due to PO4 (cation-exchange). When P-site was blocked, two basic proteins, lysozyme (Lys) and ribonuclease A(RNase), were not retained whereas cytochrome C(Cyt C) and lactoferrin (LF) were retained and also retention of acidic proteins became stronger as the repulsion due to P-site was eliminated. The number of the binding site B values determined from LGE also increased, which also showed reduction of repulsion forces. Conclusion: The selectivity (retention) of four basic proteins (RNase, Lys, Cyt C, LF) in HAC was different from that in ion-exchange chromatography. Moreover, it was possible to tune the selectivity by using NaCl gradient

HA is a calcium-phosphate complex Ca 10 (PO 4 ) 6 (OH) 2 , composed of two sites, C-site due to Ca and P-site due to PO 4 . Previous studies [8,[11][12][13]20] have shown the following general principles. i) P-site works as a cation-exchange ligand whereas Csite has a metal affinity as well as an anion-exchange ligand function.
*Address correspondence to this author at the Bio-Process Engineering Laboratory, Graduate School of Medicine, Biomedical Engineering Center (YUBEC), Yamaguchi University, Tokiwadai, Ube 755-8611, Japan; Fax +81-836-85-9201; E-mail: shuichi@yamaguchi-u.ac.jp ii) At neutral pH regions basic proteins are retained on HAC based on the electrostatic interaction with P-site and eluted with NaCl as well as sodium phosphate. This is similar to cation exchange chromatography (CIEC).
iii) Acidic proteins are retained due to a strong interaction (metal affinity) with C-site, and are not eluted with NaCl, requiring sodium phosphate for their elution. Therefore, the elution or retention behavior is different from that for anion-exchange chromatography (AIEC).
In this study, we performed linear gradient elution (LGE) of proteins on ceramic HAC columns at various conditions such as gradient slope and mobile phase pH. Several mobile phase (buffer) solution strategies and solutes were evaluated in order to probe the relative contributions of two adsorption sites of hydroxyapatite (HA) particles, C-site due to Ca (metal affinity) and P-site due to PO 4 (cation-exchange). Binding site value and elution salt concentrations determined from LGE experiments based on our method [5,25] were examined in order to elucidate the retention (binding) mechanism of HAC in comparison with the data obtained with ion-exchange chromatography (IEC).

Materials
Model protein samples used in this study are shown in Table 1. The pI values were cited from the literature [29,30]. Acidic, neutral and basic proteins of different molecular weights were chosen. Bovine serum albumin (BSA) and lysozyme (Lys) are most widely used acidic and basic model proteins, respectively. β-lactoglobulin (LG) is a good probe for the retention study as it has two genetic variants A and B forms, and the surface charge distribution is very asymmetrical [31,32]. Other reagents were of analytical grade.

Chromatography Experiments
Ceramic hydroxyapatite particles (CHT Type I, nominal particle diameter = 40 μm, Bio-Rad, USA) were self-packed into a stainless steel empty column (column inner diameter d c = 4 mm, column length Z = 100 mm, the total column bed volume V t = 1.257 mL, the column void fraction = 0.33). The void fraction ε was determined from the retention volume of high molecular weight dextran (molecular weight >2,000,000) solution pulse response curves. For comparison ion-exchange monolithic disks CIM SO3 (BIA Separations, Ajdovščina, Slovenia) were employed. They are methacrylate-based monolithic (continuous-phase) disks, in which strong cation exchange ligands (SO 3 ) were immobilized onto flow-through pores (ca. 1.5 μm) [33]. The disks (d c = 12 mm, Z = 3 mm, V t = 0.34 mL, void fraction = 0.62) were placed into the manufacturer's holder and used as a column.
Linear gradient elution (LGE) experiments were performed on a fully automated liquid chromatography system (AKTA explorer or purifier, GE Healthcare, Uppsala, Sweden) at 298±1K. Table 2 summarizes the buffer systems employed in this study. The column was equilibrated with the buffer solution until the pH and the conductivity values at the column outlet become the same as those of the buffer solution (20 to 50 column volumes). After loading the sample, the salt concentration of the mobile phase was changed with time linearly. The column was washed with the buffer solution of high salt concentrations after the sample was eluted. The column was also washed with 0.1 M NaOH once per day. The sample volume was 0.1 -0.5 mL and the sample concentration was 5-10 mg/mL. The volumetric flow-rate F was 1.0 mL/min for the CHT column and 2.0 mL/min for the CIM disks unless otherwise noted. The salt concentration at the peak retention time (or volume) I R was measured for different gradient slopes g and plotted against normalized gradient, GH. For the precise determination of I R, the dead volume between the UV detector cell and the conductivity cell was measured. The calibration curve (conductivity vs. salt concentration) was prepared for each buffer solution.
The system dead volume was reduced as much as possible. The experiments were repeated twice and the average I R value was used.
From the GH-I R curve, the number of binding sites (charges) involved in electrostatic interaction, B value, was determined according to the following equation [5,25,31,32,34]. The B value can provide valuable information on the interaction between the chromatography ligand and the solute [5,25,31,32,34]. A is the parameter including the equilibrium coefficient, the binding site and the ion-exchange capacity [5,25,31,32,34].
H is the volumetric phase ratio = (1 -ε)/ε and ε is the column void fraction = V o /V t , and gradient slope, g = (I f -I 0 )/V g , is the gradient slope of the salt (I f : final salt concentration, I 0 :initial salt concentration, V g : gradient volume). V t is the total column bed volume, V o is the column void (interstitial space) volume.

Retention Behavior as a Function of pI with Phosphate Buffers when the Two Sites are Active
Peak salt concentration I R values were measured from linear gradient elution (LGE) curves of model proteins with commonly employed phosphate buffer solutions (mobile phase A in Table 2). Relationships between I R and pI are shown in Fig. (1). It has been reported that the retention increases with pI at pH 6.8 [8]. The I R values except for βlactoglobulin, lactoferrin and cytochrome C show such a trend as shown in Fig. (1b). The I R decreased with increasing mobile phase pH, which was similar to our previous finding as shown in (Fig. 1) [25].
As already mentioned it has been reported that basic proteins are retained on P-site, which basically works as a cation-exchange ligand. Therefore, basic proteins can be eluted with NaCl as well as sodium phosphate as a gradient (modulator) component. On the other hand, acidic proteins are said to be bound to C-site based on metal affinity or coordinate binding. Acidic proteins cannot be eluted with NaCl because of its strong (affinity-based) interaction. In addition, with increasing pH the dissociation at C-site decreases and the repulsion by P-site increases.

Retention Behavior when P-site was Blocked or when NaCl was used as a Gradient Modulator
In order to examine the contribution of the two sites to the retention, LGE experiments were carried out with two different mobile phases (mobile phase B and C, Table 2). Mobile phase B uses NaCl as a gradient modulator (component) whereas Mobile phase C can suppress the P-site function by the addition of Ca ions [13]. The results with four proteins (BSA, lysozyme, lactoferrin, β-lactoglobulin) are shown in Fig. (2). BSA and lysozyme (Lys) are typical acidic and basic proteins, respectively (Fig. 1a).

Lysozyme was retained both by Mobile phase A and B
although it was not retained with Mobile phase C containing 30-fold higher Ca concentration than the former mobile phases. This is understandable as basic proteins are mainly retained on the P-site based cation exchange reaction that is blocked by Ca ions. Elution of basic proteins is possible with phosphate (Mobile phase A) or NaCl (Mobile phase B). When P-site was blocked (Mobile phase C), the retention was no longer possible.
BSA was not eluted with Mobile phase B (NaCl) due to the strong C-site based interaction. When P-site was blocked, the retention volume increased (I R increased). This is because the repulsion due to P-site was suppressed. Similar results were obtained with ovalbumin (OA).
The I R values for cytochrome C (Cyt C) and especially for lactoferrin (LF) were much higher than those for other basic proteins. Strong retention of LF was reported in the literature [35]. Similar to Cyt C and LF, β-lactoglobulin (LG) was retained much more strongly compared with other acidic proteins with Mobile phase A. Ⅱ Ⅲ Ⅳ needed. It was also retained with Mobile phase C. This indicates that LF is retained on HAC with both P-and C-sites. Cytochrome C showed similar retention behavior.
The number of binding sites (effective charges) B was determined from GH-I R curves for the three mobile phases at pH 6.8. The results are shown in (Fig. 3) along with the I R values. The B values increased for acidic proteins when Mobile phase C was employed. For two basic proteins (Lys, RNaseA), the B values decreased slightly with Mobile phase B. For another basic protein, LF the B value decreased markedly with Mobile phase B. However, the B value did not drop for Mobile phase C. Similarly, for Cyt C, the B value for Mobile C did not decrease.

Retention Behavior when C-site was Blocked
A buffer solution containing borate (1M) was proposed by Gagnon [36] to suppress C-site function although this solution was rather difficult to prepare compared with Mobile phase C. With this C-site blocking buffer (Mobile phase D, Table 2) acidic proteins (BSA, OA) and neutral protein (Mb) were not retained whereas a basic protein lysozyme was retained. However, β-lactoglobulin (LG) was retained and eluted by both NaCl and sodium phosphate gradient.
LG may have also dual interactions with C-and P-sites.

DISCUSSION
The present study has shown generally accepted retention behavior of several model proteins in HAC. Namely, basic proteins were retained due to P-site based cation exchange reaction and acidic proteins were strongly retained due to Csite. The retention became weaker with increasing mobile phase pH. However, several proteins showed different retention behavior.
Repulsion due to P-site for negatively charged biomolecules were already pointed out. For example, DNA binding is favored at high salt concentrations, which can suppress the repulsion due to P-site [20]. By suppressing the P-site function with Mobile phase C, it was possible to elucidate the binding (retention) behavior of proteins.
Four possible binding mechanisms have been proposed [19]; I) Strong affinity due to C-site, II) dual interactions with both C-and P-sites, III) C-site interaction and IV) P-site interaction. The data shown in (Fig. 1b) were tentatively classified according to these four categories.
Cytochrome C (Cyt C) and lactoferrin (LF) were retained even with Mobile phase C. It is therefore likely that both C-and P-site interactions are involved for these two proteins (category II in Fig. 1b).  Absorbance at_280 nm β-lactoglobulin (LG) was another strongly retained protein even with phosphate buffer solutions (Mobile phase A).
LG was known to be retained near the pI in both anion-and cation-exchange chromatography [31,32] due to its asymmetrical charge distribution. Arakawa and Timasheff [37] also reported this unique surface charge distribution of LG based on the solubility data. The contribution of P-site in addition to the strong interaction to the C-site might also explain the strong retention of LG.
Compared with ion-exchange chromatography (IEC), I R values are relatively low for HAC when phosphate buffer solutions are employed. At pH 6.8 I R values are below 0.1 M for most proteins. Selectivity and consequently the elution order in HAC are different from those in IEC. Moreover, in some cases as shown in (Fig. 4) the selectivity for basic proteins can be improved by using NaCl gradient instead of sodium phosphate gradient as pointed out in the literature [13][14][15][16][20][21][22].
Retention of proteins in (pure) electrostatic interactionbased chromatography (EIC) commonly called ion-exchange chromatography (IEC) is rather easy to understand. Both I R and B values increase with increasing mobile phase pH above the pI of a protein for anion-exchange chromatography (AIEC) [25,31,32]. Similarly, both I R and B values increase with decreasing mobile phase pH below the pI of a protein for cation-exchange chromatography (CIEC). Elution order of proteins basically does not change when a different IEC column is used.


extfieldsofstudy: Chemistry Medicine

text:
Simultaneous Determination of Aesculin, Aesculetin, Fraxetin, Fraxin and Polydatin in Beagle Dog Plasma by UPLC-ESI-MS/MS and Its Application in a Pharmacokinetic Study after Oral Administration Extracts of Ledum palustre L.

A rapid, simple and sensitive ultra-performance liquid chromatography-electrospray-ionization-tandem mass spectrometry (UPLC-ESI-MS/MS) method was developed and validated for the simultaneous determination of aesculin, aesculetin, fraxetin, fraxin and polydatin in beagle dog plasma for the first time. Plasma samples were pretreated by protein precipitation with methanol. Chromatographic separation was performed on an Acquity UPLC HSS T3 C18 column (2.1 mm × 100 mm, 1.8 μm) with gradient elution at a flow rate of 0.4 mL/min, using a mobile phase consisting of 0.1% formic acid (A) and acetonitrile (B). The analytes and IS were detected by multiple reaction monitoring (MRM) via negative ion mode with ion transitions of m/z 339.1–m/z 176.8 for aesculin, m/z 176.8–m/z 88.9 for aesculetin, m/z 206.8–m/z 192.1 for fraxetin, m/z 369.1–m/z 206.9 for fraxin, m/z 389.1–m/z 227.0 for polydatin and m/z 415.2–m/z 295.1 for puerarin. This method was validated according to the FDA guidelines and the results met the requirements of analysis. The calibration curves of analytes were linear with correlation coefficients more than 0.9980. The intra- and inter-day precisions were less than 15% and the accuracy was within ±15%. The maximum plasma concentration (Cmax) of aesculin, aesculetin, fraxetin, fraxin and polydatin was 46.75 ± 7.46, 209.9 ± 57.65, 369.7 ± 48.87, 67.04 ± 12.09 and 47.14 ± 12.04 ng/mL, respectively. The time to reach the maximum plasma concentration (Tmax) was 1.32 ± 0.38 h for aesculin, 1.03 ± 0.27 h for aesculetin, 0.94 ± 0.23 h for fraxetin, 0.83 ± 0.18 h for fraxin and 1.15 ± 0.15 h for polydatin. The results indicated that the absorption of aesculin might be slow in beagle dog plasma. This method was successfully applied for pharmacokinetics in beagle dog plasma after oral administration of the extracts of Ledum palustre L. at a dosage of 0.27 g/kg.

Introduction
Ledum palustre L. is an evergreen low shrub, which belongs to the Ericaceae family and is widely distributed in northern Europe, central Europe, northern Asia, America and the northeast of China. The ointment and decoction of L. palustre have been used to treat menxenia, infertility and gastric

Optimization of UPLC-ESI-MS/MS Conditions
The mobile phase was a decisive factor in the chromatographic analysis, which could affect the separation and ionization. Different solvent combinations, including water-acetonitrile and water-methanol with formic acid, were required to be tested for perfect peak shape, high response and stable retention time. It was found that the water-acetonitrile system possessed good separation efficiency with lower noise, and the retention time of the analytes was more stable than the water-methanol system. Higher signal response and sharp peaks of analytes were observed after addition of 0.1% formic acid into the mobile phase, and the sensitivity of detecting was improved. Therefore, 0.1% formic acid-acetonitrile was employed as the final mobile phase with gradient elution. The total test time was 5 min.
The analytes and IS were determined in both negative and positive ion mode, the signal sensitivity and reproducibility were better in negative ionization mode than in positive mode. Therefore, the negative ion mode was selected for the quantification. On the basis of the full-scan spectra shown in Figure 1, the deprotonated molecular ions [M − H] − were selected as parent ions because of their high response. The declustering potential (DP) and collision energy (CE) of analytes and IS were optimized in order to obtain product ion. In the process of optimization, the most abundant and stable fragment was m/z 176. 8 Table 1. and IS were optimized in order to obtain product ion. In the process of optimization, the most abundant and stable fragment was m/z 176. 8

Extraction Method
The extraction reagents play an important role in extraction recovery and matrix effect. To select the appropriate extraction method, the protein precipitation (methanol and acetonitrile) and liquidliquid extraction (LLE) (dichloromethane, ethyl acetate, acetone, mixture of ethyl acetate-acetone at different ratios) methods were evaluated by treatment of biological samples. The liquid-liquid extraction method was abandoned due to its low extraction recovery and high matrix interferences. Protein precipitation with acetonitrile has a high extraction recovery; however, strong interferences

Extraction Method
The extraction reagents play an important role in extraction recovery and matrix effect. To select the appropriate extraction method, the protein precipitation (methanol and acetonitrile) and liquid-liquid extraction (LLE) (dichloromethane, ethyl acetate, acetone, mixture of ethyl acetate-acetone at different ratios) methods were evaluated by treatment of biological samples. The liquid-liquid extraction method was abandoned due to its low extraction recovery and high matrix interferences. Protein precipitation with acetonitrile has a high extraction recovery; however, strong interferences were observed. The recovery of analytes and IS were higher and more stable with methanol, with no obvious interferences. Therefore, methanol was chosen as the extraction solvent.

Selectivity and Specificity
The six beagle dogs' blank plasma, spiked plasma samples with the analytes and IS at lowest limit of quantification (LLOQ) level, spiked plasma samples with the analytes and IS at quality control middle (QCM) level, and plasma samples 1 h after oral administration of the extract of L. palustre are displayed in Figure 2. The retention time of aesculin, aesculetin, fraxetin, fraxin, polydatin and IS were 2.96, 3.27, 3.34, 3.10, 3.48 and 3.07 min, respectively. No apparent interferences were observed from the endogenous substance during the retention time of the analytes and IS.

Linearity and Sensitivity
The linearity status of analytes such as typical regression equations, correlation coefficients and linear ranges are displayed in

Precision and Accuracy
The measurement results for precision and accuracy on LLOQ and QC samples (low, middle and high) are shown in Table 3. The intra-day and inter-day precisions of the analytes were lower than 16.1%, and the accuracy ranged from −5.2% to 16.5%. All values were within acceptable limits, indicating that the method described above has better repeatability and reproducibility.

Extraction Recovery and Matrix Effect
The extraction recovery of the five analytes at low, middle and high concentrations ranged from 80.10% to 98.54%, with RSD < 10.2%. The results indicate that the extract method of protein precipitation is suitable for the determination. The matrix effects of the analytes at the three concentrations were within the range of 95.86% to 106.28%, with RSD < 13.2%, which revealed no obvious ion suppression or enhancement for the analytes in beagle dog plasma. The extraction recovery and matrix effect of IS were 96.81% and 98.74%, respectively. The data obtained are presented in Table 4.

Stability
The test data for the stability of the analytes at three concentrations under various conditions are presented in Table 5. The results indicate that all samples were stable under the conditions mentioned above, and that the stability of samples had not been affected during the analysis process.

Application to a Pharmacokinetic Study
The validated UPLC-ESI-MS/MS method was successfully applied for the simultaneous determination of aesculin, aesculetin, fraxetin, fraxin and polydatin in beagle dog plasma after oral administration of the extract of L. palustre at a dosage of 0.27 g/kg. The pharmacokinetic parameters analyzed are presented in Table 6, and the mean plasma concentration-time profiles are illustrated in Figure 3. The LLOQ of aesculin was 0.50 ng/mL, 1.09 ng/mL for aesculetin, 2.15 ng/mL for fraxetin, 1.25 ng/mL for fraxin and 0.50 ng/mL for polydatin. The concentration of Cmax-to-LLOQ ratio of the analytes was greater than 20. The results indicated that the UPLC-MS/MS method has adequate sensitivity to measure the five components in beagle dog plasma.   The content of the five analytes was determined by HPLC in this paper. We can conclude that the content of aesculin was higher than aesculetin, and that fraxin was higher than fraxetin in the extract of L. palustre. However, the value of the last measurable plasma concentrations (AUC0 → t) of aesculetin and fraxetin exceed those of aesculin and fraxin in Table 6. The maximum plasma concentrations (Cmax) of aesculin, aesculetin, fraxetin, fraxin and polydatin were 46.75 ± 7.46, 209.9 ± 27.65, 369.7 ± 48.87, 67.04 ± 12.09 and 47.14 ± 12.04 ng/mL, respectively; it can be seen that the Cmax of fraxetin was higher than others. The phenomena above may be due to the fact that aesculin and fraxin can be transformed into aesculetin and fraxetin, respectively, in beagle dog plasma. According to the chemical structures of the analytes in Figure 1, aesculin and fraxin are the glucosides of aesculetin and fraxetin. Previous studies have revealed that aesculin can be completely metabolized to aesculetin under the action of gut bacteria, because it degrades the glycoside of aesculin to facilitate its absorption in the intestine [23]. Similar experimental results have also appeared in other reports about the metabolic fate of fraxinin, indicating that fraxin can be extensively metabolized to fraxetin in vivo, and that intestinal microflora play an important role in this process [24]. The time to reach the maximum plasma concentration (T max ) was 1.32 ± 0.38 h for aesculin, 1.03 ± 0.27 h for aesculetin, 0.94 ± 0.23 h for fraxetin, 0.83 ± 0.18 h for fraxin and 1.15 ± 0.15 h for polydatin. These results indicate that the absorption of aesculin might be slow in beagle dog plasma. The elimination half-life (T 1/2 ) was 3.43 ± 0.47 h for aesculin, 4.25 ± 0.18 h for aesculetin, 3.76 ± 0.35 h for fraxetin, 2.99 ± 0.29 h for fraxin and 3.36 ± 0.31 h for polydatin. These results indicate that the elimination of fraxetin may be slow because of the longer T 1/2 than the other four components.
The content of the five analytes was determined by HPLC in this paper. We can conclude that the content of aesculin was higher than aesculetin, and that fraxin was higher than fraxetin in the extract of L. palustre. However, the value of the last measurable plasma concentrations (AUC 0→t ) of aesculetin and fraxetin exceed those of aesculin and fraxin in Table 6. The maximum plasma concentrations (C max ) of aesculin, aesculetin, fraxetin, fraxin and polydatin were 46.75 ± 7.46, 209.9 ± 27.65, 369.7 ± 48.87, 67.04 ± 12.09 and 47.14 ± 12.04 ng/mL, respectively; it can be seen that the C max of fraxetin was higher than others. The phenomena above may be due to the fact that aesculin and fraxin can be transformed into aesculetin and fraxetin, respectively, in beagle dog plasma. According to the chemical structures of the analytes in Figure 1, aesculin and fraxin are the glucosides of aesculetin and fraxetin. Previous studies have revealed that aesculin can be completely metabolized to aesculetin under the action of gut bacteria, because it degrades the glycoside of aesculin to facilitate its absorption in the intestine [23]. Similar experimental results have also appeared in other reports about the metabolic fate of fraxinin, indicating that fraxin can be extensively metabolized to fraxetin in vivo, and that intestinal microflora play an important role in this process [24].

Instruments and Chromatographic Conditions
A Waters Acquity UPLC system (Waters, Milford, MA, USA) was used for liquid chromatography analysis. Chromatographic separation was performed on an Acquity UPLC HSS T3 C18 column (2.1 mm × 100 mm, 1.8 µm) with a mobile phase composed of 0.1% formic acid (A) and acetonitrile (B). The gradient elution mode at a flow rate of 0.4 mL/min was as follows: 0 min at 2% B, 0-3 min at 2%-80%, 3-5 min at 80%-2%. The column temperature was maintained at 25 • C. Samples were analyzed by injecting 5 µL into the instrument, with a total run time of 5 min.
An  Figure 1. Mass spectrometric parameters were optimized for better ionization of analytes. The optimal MS parameters were as follows: source temperature at 450 • C, capillary voltage at 3200 V, Nitrogen was used as desolvation and cone gas with a flow rate of 700 L/h and 150 L/h, cone voltage at 30 V.

Preparation of L. palustre Extract
The L. palustre was extracted by refluxing with ethanol three times, each time for 1 h, then filtrated with gauze. The filtrate was combined and concentrated up to dryness. Based on measurement, we could see that the extraction rate was 23%. According to the conversion relationships between humans and animals, the oral dosage for dogs was determined to be 0.27 g/kg. The contents of aesculin, aesculetin, fraxetin, fraxin and polydatin in L. palustre were determined by HPLC (Waters). The content of aesculin was 5.9 mg/g, 0.9 mg/g for aesculetin, 2.1 mg/g for fraxetin, 2.3 mg/g for fraxin and 7.4 mg/g for polydatin in extracts of L. palustre, respectively. The extracts were fed to beagle dogs in the form of capsules.

Preparation of Calibration Standards and Quality Control (QC) Samples
A mixed stock solution was prepared in methanol and the concentration of the analytes was 97.50 ng/mL for aesculin, 209.2 ng/mL for aesculetin, 412.8 ng/mL for fraxetin, 241.0 ng/mL for fraxin,

Sample Preparation
Protein precipitate method was used for treatment of the biological plasma samples. 100 µL of IS and 300 µL of methanol were added to a 100 µL plasma sample in a 10 mL glass tube. The tube was vortexed (Qingpu Huxi instrument Factory, Shanghai, China) for 3 min and centrifuged (Anhui USTC Zonkia Scientific Tnstruments Co., Ltd., Anhui, China) at 3500 rpm for 10 min. Supernatant was transferred to a new glass tube and evaporated to dryness under nitrogen stream at 40 • C. The residue was reconstituted with 100 µL of initial mobile phase, then vortexed for 3 min and centrifuged at 12,000 rpm for 10 min. A 5 µL aliquot was injected into the UPLC-ESI-MS/MS system for analysis.

Method Validation
Method validation was implemented in beagle dog plasma according to the US Food and Drug Administration Bio-analytical Method Validation Guide [25]. The test projects include selectivity, linearity, sensitivity, precision, accuracy, matrix effect, recovery, stability.

Selectivity and Specificity
The specificity of the method was estimated by comparing the chromatograms of the blank plasma from six beagle dogs, spiked plasma samples with the analytes and IS at LLOQ level, spiked plasma samples with the analytes and IS at QCM level, and plasma samples from beagle dog 1 h after oral administration of the extract of L. palustre. The purpose of specificity research was to detect whether endogenous interference from other substances has an effect on retention time of analytes and IS.

Linearity and Sensitivity
The linearity of the five compounds was evaluated via construct regression equation. Calibration curve was established by plotting peak area ratios of the analytes to IS versus plasma at eight concentration levels in the linear range of analytes, using the least-square linear regression with weighting factor 1/X 2 . The regression equation was Y = aX + b.
The lowest limit of quantification (LLOQ) was defined as the lowest concentration in the calibration curve at which the signal-to-noise (S/N) ratio was about 10. The precision and accuracy of LLOQ should be within acceptable limits. The precision was required to be less than 20% and the accuracy was required to be less than ±20%.

Precision and Accuracy
The precision and accuracy were evaluated by analyzing six replicates of QC samples at three concentration levels. The intra-day and inter-day precisions were determined on the same day and on three consecutive days. Precision and accuracy were represented by RSD (%) and RE (%), respectively. Furthermore, the values of precision were required not to exceed 15%, and the accuracy should be better than ±15%.

Extraction Recovery and Matrix Effect
The extraction recovery and matrix effect were determined in six replicates at three QC concentrations. Extraction recovery of analytes was evaluated by comparing the peak areas of the analytes from the extracted QC samples with those samples spiked after extraction of blank plasma at equivalent concentrations. The matrix effect of the analytes was assessed by comparing the peak areas of the analytes obtained from the plasma samples with those extracted from standard solutions at the same concentration. The extraction recovery and matrix effect of IS was also determined in a similar way at a single concentration of 260.0 ng/mL.

Stability
The stability of the analytes in beagle dog plasma was measured by analyzing six replicates of QC samples at three concentration levels in various storage conditions, including short-term stability (kept at room temperature for 24 h), long-term stability (kept at −20 • C for 2 weeks), three freeze/thaw cycles stability (repeated freeze/thaw for three cycles from −20 • C to room temperature) and post-preparation stability (processed samples kept under the auto-sample conditions for 24 h). Samples were considered stable if the deviation was within 15% compared with a freshly prepared solution.  ) with a non-compartmental model was used to acquire pharmacokinetic parameters. The values of maximum plasma concentration (C max ) and time to reach the maximum plasma concentration (T max ) were observed from the mean plasma concentration-time profiles. The terminal elimination rate constant (K e ) was measured by linear regression of the terminal log-linear part of the mean plasma concentration-time curve, and the T 1/2 = 0.693/K e was used to measure elimination half-life (T 1/2 ). The area under the plasma concentration-time curve at the time from 0 to the last measurable plasma concentration (AUC 0→t ) was estimated by using the linear trapezoidal rule. The formula of AUC 0→∞ = AUC 0→t + C t /K e can be used for calculating the area under the plasma concentration-time curve at any time from 0 to infinity (AUC 0→∞ ). All analysis results were expressed as mean ± SD.

Conclusions
In this paper, a rapid, sensitive and selective UPLC-ESI-MS/MS method was successfully applied for simultaneous determination of aesculin, aesculetin, fraxetin, fraxin and polydatin in beagle dog plasma. This method provided adequate recovery and matrix effect with precision, accuracy and lower detection limits. All test results are within acceptable limits with a total run time of 4.5 min. This is the first report about the pharmacokinetics of aesculin, aesculetin, fraxetin, fraxin and polydatin in beagle dog plasma. Therefore, the information obtained provides a reference for further clinical application of L. palustre.