Jhimli Mitra
email: jhimlimitra@yahoo.com

Robert Martí

Arnau Oliver
email: aoliver@eia.udg.edu

Xavier Lladó
email: llado@eia.udg.edu

Soumya Ghose
email: soumyaghose@gmail.com

Joan C Vilanova
email: kvilanova@comg.cat

Fabrice Meriaudeau
email: fabrice.meriaudeau@u-bourgogne.fr

B-splines Coupled with Quadrature Texture to Register Prostate Multimodal Images

Keywords: Prostate multimodal registration, Bspline, quadrature filter, texture energy, normalized mutual information

Purpose: Needle biopsy of the prostate is guided by Transrectal Ultrasound (TRUS) imaging. The TRUS images do not provide proper spatial localization of malignant tissues due to the poor sensitivity of TRUS to visualize early malignancy. Magnetic Resonance Imaging (MRI) has been shown to be sensitive for detection of early stage malignancy and therefore, a novel 2D deformable registration method that overlays pre-biopsy MRI onto TRUS images has been proposed. Method: The registration method involves B-spline deformations with Normalized Mutual Information (NMI) as similarity measure computed from the texture images obtained from the amplitude responses of the directional quadrature filter pairs. Registration accuracy of the proposed method is evaluated by computing the Dice Similarity coefficient (DSC) and 95% Hausdorff Distance (HD) values for 20 patient datasets and Target Registration Error (TRE) for 18 patients only where homologous structures are visible in both the TRUS and transformed MR images. Results: The proposed method and B-splines using NMI computed from intensities provide average TRE values of 2.64 ± 1.37 mm and 4.43 ± 2.77 mm respectively. Our

Introduction

Transrectal Ultrasound (TRUS) is routinely used for interventional needle biopsy of the prostate to detect prostate cancer. However, ultrasound images only provide very limited information about the location of prostate cancer. Cancer tissues are either hypoechoic or isoechoic in TRUS. The isoechoic cancer tissues account for about 24% -42% [START_REF] Carroll | Transrectal ultrasound guided prostate biopsy[END_REF]. As reported in [START_REF] Arcangeli | Prostate-specific antigen as a screening test for prostate cancer-The United States experience[END_REF] the probability to diagnose prostate cancer from TRUS biopsy alone is about 70% -80%. Therefore, the multicore biopsy technique used for needle biopsy oftens fails to extract positive samples and the number of rebiopsies thereby increases. It has been demonstrated in [START_REF] Vilanova | Usefulness of prebiopsy multi-functional and morphologic MRI combined with the free-to-total PSA ratio in the detection of prostate cancer[END_REF] that detection of prostate cancer using Magnetic Resonance Imaging (MRI) has a negative predictive value of 95% or greater and also the accuracy of MRI to diagnose prostate cancer is 95%. Therefore, MRI may serve as a triage test for men deemed to be at risk of prostate cancer. In other words, if MRI has a negative predictive value ≥ 95% for significant cancer then it might greatly reduce the number requiring biopsy, while at the same time increasing the yield of useful information for those who are sent for biopsy. Therefore, fusion of pre-biopsy MR images onto interventional TRUS images might increase the overall biopsy accuracy [START_REF] Kaplan | Real time MRI-ultrasound image guided stereotactic prostate biopsy[END_REF], [START_REF] Xu | Real-time MRI-TRUS fusion for guidance of targeted prostate biopsies[END_REF], [START_REF] Singh | Initial clinical experience with real-time transrectal ultrasonographymagnetic resonance imaging fusion-guided prostate biopsy[END_REF]. Fei et al. [START_REF] Fei | Slice-to-volume registration and its potential application to interventional MRI-guided radio-frequency thermal ablation of prostate cancer[END_REF] registered interventional Magnetic Resonance Imaging (iMRI) slices to previously obtained high resolution MRI volume of prostate images for iMRI guided radio frequency thermal ablation of prostate cancer. A validation and integration for 3D TRUS guided robotic surgery for prostate brachytherapy was done by Wei et al. [START_REF] Wei | Robot-assisted 3D-TRUS guided prostate brachytherapy: System integration and validation[END_REF]. Hummel et al. [START_REF] Hummel | 2D/3D registration of endoscopic ultrasound to CT volume data[END_REF] proposed a rigid registration for fusion of 2D intraoperative US with preoperative CT volume based on fiducials.

The prostate of the same patient may undergo deformations during the TRUS or MR imaging procedures due to full bladder or rectum, altered patient positions on the couch and inflation of the endorectal balloon inside the rectum containing the MR coil. Intended to cope with these deformations, nonrigid registration methods are applied for prostate multimodal fusion [START_REF] Mizowaki | Towards integrating functional imaging in the treatment of prostate cancer with radiation: The registration of the MR spectroscopy imaging to ultrasound/CT images and its implementation in treatment planning[END_REF], [START_REF] Reynier | MRI/TRUS data fusion for prostate brachytherapy. preliminary results[END_REF], [START_REF] Du Bois D'aische | Improved non-rigid registration of prostate MRI[END_REF], [START_REF] Narayanan | MRI-ultrasound registration for targeted prostate biopsy[END_REF], [START_REF] Chen | 3D meshless prostate segmentation and registration in image guided radiotherapy[END_REF], [START_REF] Baumann | Prostate biopsy assistance system with gland deformation estimation for enhanced precision[END_REF]. Spline based deformations have been commonly used for prostate registration. Lu et al. [START_REF] Lu | A statistical volumetric model for characterization and visualization of prostate cancer[END_REF] used Thin-plate Splines (TPS) warping to generate statistical volumetric model of prostate for localization of prostate cancer. However, the Target Registration Error (TRE) value of 295.66 pixels seemed to be too large to be used in clinical interventions. Prostate MR volumes were warped using TPS by Fei et al. [START_REF] Fei | A comparative study of warping and rigid body registration for the prostate and pelvic MR volumes[END_REF] where the prostate centroid displacement was observed to be 0.6 mm. Vishwanath et al. [START_REF] Vishwanath | Integrating structural and functional imaging for computer assisted detection of prostate cancer on multi-protocol In Vivo 3 Tesla MRI[END_REF] registered prostate histological slices and MR slices to detect prostate cancer using B-splines. Since the aim was to detect cancer, quantitative values related to registration accuracy were not presented. A recent work by Xiao et al. [START_REF] Xiao | A structuralfunctional MRI-based disease atlas:application to computer-aided-diagnosis of prostate cancer[END_REF] proposed to build a spatial disease atlas of the prostate using both B-splines and TPS. However, only qualitative results were presented. Oguro et al. [START_REF] Oguro | MRI signal intensity based B-spline nonrigid registration for pre-and intraoperative imaging during prostate brachytherapy[END_REF] registered pre-and intraoperative MR images for prostate brachytherapy. Dice Similarity Coefficient (DSC) value for the total gland was reported as 0.91 and the fiducial registration error was 2.3 ± 1.8 mm respectively. TRUS and MR multimodal registration for interventional biopsy was attempted by Mitra et al. [START_REF] Mitra | A thin-plate spline based multimodal prostate registration with optimal correspondences[END_REF]. The method was based on TPS with automatic point correspondences generated with an average DSC of 0.97 ± 0.01 for 4 patient datasets.

In this work a method to register TRUS and MR prostate 2D images with B-spline free-form deformations and quadrature filter pairs is presented. The Bspline registration uses uniform grids over the MR image and NMI is used as a similarity measure [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF], [START_REF] Ino | A data distributed parallel algorithm for nonrigid image registration[END_REF]. However, the novelty of the proposal is to use directional quadrature filter pairs to transform both the MR and TRUS images into texture images obtained from the amplitude response of the filter pairs and use these transformed images for NMI computation. A similar method has been used by Jarc et al. [START_REF] Jarc | Texture features for affine registration of thermal (FLIR) and visible images[END_REF] employing Law's texture to compute the MI in order to register far-infrared and visible spectrum gray-scale images. Francois et al. [START_REF] Francois | Robust statistical registration of 3D ultrasound images using texture information[END_REF] used texture-based statistical measures to register carotid ultrasound volumes where the texture information was given by spatial Gabor filters. The remaining paper is organized with the details of the method in section 2. Results and discussions follow in section 3 with the conclusions in section 4.

Materials and Methods

The Proposed Method and Algorithm

The proposed method uses B-spline deformations to register 2D MR slice corresponding to the 2D TRUS slice. This work is a feasibility study (improving on registration accuracy) on the type deformable registration algorithm that may be used during prostate biopsy for accurate extraction of biopsy samples. Therefore, the corresponding 2D axial MR slice of the respective axial TRUS slice is identified by an expert that closely match each other. In the clinical context, the idea would be to track the z-coordinate of the TRUS axial 2D slice with an electromagnetic tracking device attached to the TRUS probe during biopsy and find an automatic slice correspondence with a reconstructed 3D TRUS volume and finally register to the 3D MR volume [START_REF] Xu | Real-time MRI-TRUS fusion for guidance of targeted prostate biopsies[END_REF]. In our work, once the slice correspondences are established, the prostate in each of the images is then manually segmented. However, we plan to use some automatic segmentation methods proposed by Shen et al. [START_REF] Shen | Segmentation of prostate boundaries from ultrasound images using statistical shape model[END_REF], Cosio [START_REF] Cosío | Automatic initialization of an Active Shape Model of the prostate[END_REF] and Ghose et al. [START_REF] Ghose | Texture guided Active Appearance Model propagation for prostate segmentation[END_REF] for automatic real-time fusion of multimodal prostate images. The optimal B-spline deformations are obtained by maximizing the NMI as similarity measure. The NMI of the fixed TRUS and the moving MR images is computed from the image texture energies obtained from the amplitude response of the directional quadrature filters (section 2.3). The advantage of computing NMI by this method is that the gray level differences in the TRUS and MR modalities are minimized, enhancing the inherent texture information of the prostate. A schematic diagram in Fig. 1 shows the work flow of the proposed algorithm.

Registration using B-splines

A common technique to represent a free-form deformation is to employ spline functions as B-splines [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF]. B-splines consist of a set of control points that can be locally controlled on the image domain. Although, the original B-splines equation is formulated for 3Ddeformable objects [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF], [START_REF] Ino | A data distributed parallel algorithm for nonrigid image registration[END_REF], the spline functions in this paper are represented for 2D images.

Let Ω = {(x, y)|0 ≤ x < X, 0 ≤ y < Y } represent the image domain. The transformation between the moving and fixed images is given by T:(x, y) → (x ′ , y ′ ), where any point (x, y) of the moving image is mapped onto its corresponding point (x ′ , y ′ ) on the fixed image. Given a mesh of control points on the moving image with a control point defined as φ i,j with uniform spacing of δ mm, the nonrigid transformation T is defined by B-spline functions as

T(x, y) = 3 l=0 3 m=0 B l (u)B m (v)φ i+l,j+m (1)
where i = ⌊x/δ⌋ -1, j = ⌊y/δ⌋ -1, u = x/δ -⌊x/δ⌋ and v = y/δ -⌊y/δ⌋. ⌊.⌋ is the floor function and B l represents the l th basis function of the cubic B-spline functions such that

B 0 (u) = (1 -u 3 )/6 B 1 (u) = (3u 3 -6u 2 + 4)/6 B 2 (u) = (-3u 3 + 3u 2 + 3u + 1)/6 (2) B 3 (u) = u 3 /6.
The B-spline free-form deformations are locally controlled because the deformation at any point (x, y) is controlled by its neighboring 4 × 4 control points. Bsplines provide a wide range of deformations by organizing the mesh of control points and the images in a hierarchy [START_REF] Forsey | Hierarchcal B-spline refinement[END_REF], i.e. the distance between the control points decrease introducing more control points while the images move from coarser to finer levels using a multiresolution Gaussian pyramid. The B-spline control points grid refinement is done using the standard splitting matrix [START_REF] Yamaguchi | Curves and Surfaces in Computer Aided Geometric Design[END_REF].

The similarity measure used for B-splines deformation is NMI (described in details in section 2.4) between the moving (M ) and the fixed (F ) images and the optimization is solved using a quasi-Newton optimization method as "Limited Memory Broyden-Fletcher-Glodfarb-Shanno" (L-BFGS) algorithm [START_REF] Liu | On the limited memory method for large scale optimization[END_REF].

Quadrature Filters

Band-pass quadrature filters have been used in computer vision to access multi-scale image information like local-phase, energy, angular frequency etc. Central to the theory of quadrature filters in the analytical domain is the Hilbert transform [START_REF] Boukerroui | On the choice of band-pass quadrature filters[END_REF]. The analytical signal of a 1D real signal f (x) is given by

f A (x) = f (x) -if H (x); (3)
where i = √ -1 and f H (x) is the Hilbert transform of f (x) defined by:

f H (x) = 1 π +∞ -∞ f (τ ) τ -x dτ (4)
⇔ F H (ω) = F (ω).i sign(ω),
where, F (ω) is the Fourier transform of f (x) and

sign(ω) = -1 ω < 0 +1 ω ≥ 0 (5)
Therefore, the analytical signal in Fourier domain is obtained from (3) and ( 5) as

F A (ω) = F (ω). [1 + sign(ω)] . (6)
To compute the local features of an image, localization of both space and frequency is required and is not possible directly from the resulting analytical signal since the Hilbert transform or analytical signal in (3) and ( 5) is defined over the entire signal. Therefore, an alternative approach is to compute the local phase or energy from the filtered version of the signal. The filter is an even, symmetric, zero-DC filter f e (x) and its odd counterpart is f o (x) which is the Hilbert transform of f e (x), hence they are in quadrature. Therefore, the analytical signal can be written as Where H(.) is the Hilbert transform and '*' is the 1-D convolution operator.

fA (x) = f e (x) * f (x) -iH(f e (x) * f (x)) (7) = (f e (x) -iH(f e (x)) * f (x) = (f e (x) -if o (x)) * f (x).
In practice, an approximation of the local amplitude or energy ( Â(x)) and phase ( φ(x)) is obtained by using band-pass quadrature even and odd filter pair, f e (x) and f o (x) respectively, where,

Â(x) = [f e (x) * f (x)] 2 + [f o (x) * f (x)] 2 (8) φ(x) = arctan {f e (x) * f (x)/f o (x) * f (x)} . (9)
A generalization of quadrature filters in 2-D is provided as a set of filters tuned to a particular orientation and are therefore called directional quadrature filters [START_REF] Jähne | Digital Image Processing[END_REF].

The magnitude of the complex filter response with images gives the power portions of the texture process contained in different spectral bands [START_REF] Aach | On texture analysis: Local energy transforms versus quadrature filters[END_REF]. Log-Gabor quadrature filters are used in this paper that are Gaussian functions on logarithmic scale. The 1-D representation of log-Gabor function in frequency domain is given by

G l (ω) = exp - ln 2 (ω/ω 0 ) 2 ln 2 (κ β ) ( 10
)
where ω 0 is the peak tuning frequency at π/3 and 0 < κ β < 1.

κ β = exp - 1 4 2 ln(2)β ,
where β is the bandwidth fixed to 2 octaves. The peak tuning frequency and bandwidth are optimized values obtained from the MATLAB toolbox provided by Andersson and Knutsson [START_REF] Andersson | Adaptive filtering[END_REF]. Quadrature filters in directions 0 • , 45 • , 90 • and 135 • are used and the individual even and odd filter responses are added to provide the magnitude of the combined filter response. The even and odd filters tuned to 0 • are shown in Fig. 2. Fig. 3 shows the texture energies of the fixed TRUS and the moving MR images when the 4 directional even-odd quadrature filter pairs are applied. It is evident that the gray level differences between the internal structures of the prostate are minimized in MR image as well as the shadow regions have disappeared in the TRUS image.

Normalized Mutual Information

The magnitude responses of the quadrature filters obtained in section 2.3 are used to compute the NMI between the moving and fixed images. The NMI is an information theoretic measure that tries to reduce the joint entropy of the images [START_REF] Studholme | An overlap invariant entropy measure of 3D medical image alignment[END_REF] and is given by

NMI = ζ similarity = H(M ) + H(F ) H(M, F ) (11)
where ζ similarity is the similarity measure for B-splines registration that is maximized in the process, H(M ) and H(F ) are the marginal entropies of the moving (M ) and fixed (F ) images respectively, and H(M, F ) is the joint entropy of the images. H(M, F ) can be written using probability theory as

H(M, F ) = - m,f p(m, f ) log [p(m, f )], (12)
where, p(m, f ) is the joint probability distribution of the images obtained from their joint histogram.

Data and Analysis

The proposed method is evaluated with datasets of 20 patients with average prostate volume of 56. Common metrics to evaluate registration accuracies are DSC [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF], TRE [START_REF] Maurer | Estimation of accuracy in localizing externally attached markers in multimodal volume head images[END_REF], [START_REF] Maurer | Registration of head volume images using implantable fidu-cial markers[END_REF] and 95% Hausdorff Distance (HD) [START_REF] Huttenlocher | Comparing images using the Hausdorff distance[END_REF]. DSC is a measure of overlap of the same labels (E) between the transformed moving image (M (E)) and the fixed image (F (E)) and is given by

DSC = 2(M (E) ∩ F (E)) M (E) + F (E) . ( 13
)
A target is an anatomical landmark in the patient's body and is normally the centroid of a lesion, tumor, gland, etc. that is not used to compute the transformation of the moving image to the fixed image. TRE is the root mean square distance of such homologous targets tp i and tq i , i = 1, 2, • • • , N on the moving and the fixed images respectively and is given by

TRE = 1 N N i=1 (T(tp i ) -tq i ) 2 (14)
where, T(.) is the transformation of the moving image. The targets used in our experiments are primarily centroids of lesions and tumors in central gland, prostatic urethra, sometimes centroids of tumor in peripheral region and the centroid of the central gland in few cases where lesions or other homologous structures are not visible in TRUS as in the corresponding MRI. The repeatability error in the localization of the targets is given as the Target Localization Error (TLE) computed from the centroids of manually selected regions from 5 independent trials by an experienced radiologist. Given a set of finite points A = {a 1 , . . . , a p } and B = {b 1 , . . . , b q }, the Hausdorff distance between the point sets is defined by

HD(A, B) = max(h(A, B), h(B, A)) (15)
where

h(A, B) = max a∈A (min b∈B a -b ) ( 16
)
and . is the Euclidean norm.

We have compared the proposed method against two spline-based registration methods; 1) B-splines registration that maximizes the NMI computed from the raw intensities of the multimodal images [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF], and 2) registration using TPS that uses contour-based automatic correspondences to solve the affine and TPS weight parameters [START_REF] Mitra | A thin-plate spline based multimodal prostate registration with optimal correspondences[END_REF].

Results and Discussions

The registration accuracies of our method are evaluated in terms of DSC and 95% HD values for 20 patients and TRE and TLE values for 18 patients i.e. only where homologous structures are visible both in TRUS and the transformed MR images. Table 1 shows the DSC, TRE, TLE and 95% HD values as obtained from the experiments for B-spline deformations with NMI computed from intensities [START_REF] Rueckert | Nonrigid registration using free-form deformations: Application to breast MR images[END_REF], our method for B-spline deformations with NMI computed from textures and TPS registration [START_REF] Mitra | A thin-plate spline based multimodal prostate registration with optimal correspondences[END_REF]. A statistically significant reduction with two-tailed t-test p = 0.02 is observed in the average TRE value for B-splines with NMI computed from texture when compared to B-splines with NMI computed from intensities with 2.64 ± 1.37 mm and 4.43 ± 2.77 mm respectively. Our proposed method shows an improvement of 1.18 times in TRE when compared with TPS registration with average TRE of 3.11 ± 2.18 mm. However, this improvement of average TRE over TPS is not statistically significant. The average TLE for all the methods are similar with 0.16±0.11 mm, 0.12±0.08 mm and 0.15 ± 0.12 mm for our method, B-splines with NMI from intensities and TPS registration respectively. It is observed from Table 1 that TLE is normally higher for patient cases where the centroid of a larger region is considered as the target.

The DSC value is a global measure of region overlap and the average DSC values for our method of B-splines with NMI computed from texture, B-splines with NMI computed from intensities and TPS registration are 0.943 ± 0.039, 0.944 ± 0.036 and 0.975 ± 0.007 respectively. 95% HD provides the contour accuracy for which the average values of 4.75 ± 3.40 mm, 4.85 ± 2.30 mm and 2.17 ± 0.85 mm are obtained for our method of B-splines with NMI from texture, B-splines with NMI from intensties and TPS respectively. We observe that TPS registration based on control points placed over the contours always has higher average DSC and lower average HD values over the proposed B-splines registration with statistical significance (two-tailed t-test) of p < 0.001 and p < 0.0001 for DSC and 95% HD measures respectively. TRE values of all patients are much lower with the proposed method when compared to the B-spline method with NMI computed from intensities except for patients 2, 7, 13 and 18 where the dark shadows near the edges of the TRUS images are misinterpreted as the black background. Higher DSC values are obtained from the proposed method when compared to B-spline with NMI computed from intensities especially for patients 6, 9, 10, 14, 15, 19 and 20 and for patients 2, 9, 10 and 20 when compared to TPS. Fig. 5 shows the qualitative B-spline registration using NMI from raw intensities, using NMI from texture and the TPS registration in the form of checkerboards for 6 patients (patient 6, 8, 10, 11, 13 and 15). As seen in Fig. 5 for patients 6 (column 1), 10 (column 3) and 15 (column 6), the checkerboards (row 6) show good region overlaps that are also evident from Table 1 for the proposed method. TRE is a reliable measure of registration accuracy than region overlap measure specially when localization of biopsy site is involved. Therefore, in spite of the less satisfactory region overlaps in Fig. 5 for the proposed method (row 6) with patients 8 (column 2), 11 (column 4) and 15 (column 6) when compared to other methods (rows 4 and 8), lower TRE values are obtained for the same patient cases with our proposed method as seen in Table 1. In Fig. 5 we observe poor region overlaps for patients 11 and 13 in columns 4 and 5 respectively. There may be two reasons for the poor registration accuracy around the prostate contour, 1) acoustic shadows in TRUS images around the rectum do not provide any texture information and are considered homogeneous with the black background, and 2) part of the contour of the moving image may lie far from the respective fixed image contour. In both these cases, the control grids placed on the moving image around the contour consider a large part of the black background for the maximization of NMI with the corresponding TRUS (textured) region and therefore the maximization process fails to reach a global maximum. However, blurring of the prostate around the bladder or rectum would not affect the contour registration accuracy since blurred regions still may contain some texture. As seen for patient 15 (row 1, column 6) that the blurring of the prostate around the rectum in TRUS does not affect the overlap accuracy. It is also observed that shadow artifacts and calcifications inside the prostate do not affect the contour registration accuracy. For instance, for patient 11 (row 1, column 4), the acoustic shadow region inside the prostate does not affect the overlap accuracy; however, the large part of the shadow on the lower-right contour region deteriorates the contour overlap accuracy.

The texture obtained from the magnitude of directional quadrature filter transforms of an image measures the power portions of the image. The MR and TRUS images have varied gray-level intensities and contrasts. Therefore, transforming the MR and TRUS images as texture energy images homogenizes the intensity variations between them and reveals the underlying prostate architectural information that is common to both the modalities. Hence, the proposed method shows better registration accuracies in terms of TRE than traditional B-splines deformation with NMI from intensities. The proposed algorithm is validated only on TRUS-MR slices from the prostate mid-gland region and the performance of the same on the base and apex regions is yet to be validated.

The algorithms are implemented in MATLAB 2009b with a machine configuration of 1.66 MHz Core2Duo processor and 2 GB memory. The average time requirements of our method of B-splines with NMI computed from texture, B-splines with NMI computed from intensities and TPS registrations are 797.72±202.59 secs, 147.25 ± 43.81 secs and 76.22 ± 29.79 secs respectively. The computation time includes the time for affine and nonrigid registrations for each of the methods. The obvious reason of B-splines with NMI computed from texture being computationally expensive is due to the use of 4 quadrature convolutions for each of the fixed and moving images at each resolution of B-splines deformations. The lowest time of 397 secs and the maximum time 1316 secs are recorded for the proposed method. However, quadrature convolution and registration of 256 × 256 × 256 3D volumes are achievable at 3.05 secs when programmed in GPU [START_REF] Eklund | Phase based volume registration using CUDA[END_REF]. Therefore, parallelization of the convolutions and our algorithm implemented on GPU would ideally reduce the execution time to less than 3 secs which is closer to the clinical requirement.

Conclusions

A method to register TRUS and MR prostate 2D images have been presented that uses B-spline deformations with a novel method of computing the NMI. The NMI as similarity measure for the registration is computed from texture energy of the images obtained from magnitude of the directional quadrature filter pair responses. Log-Gabor filters with narrow bandwidth has been used that allows to measure power portions of the signal representing texture energy in case of 2D TRUS and MR images. NMI computation involves reduction of entropy of the images. The entropy between TRUS and MR raw intensity images is typically more than the entropy of texture images due to variations in the gray-levels. Therefore, B-spline registration with NMI computed from texture images provides more average accuracy in terms of TRE than that with NMI computed from intensity images. We observed from the experiments that the average DSC and 95% HD values for TPS registration show smaller error compared to the proposed method and B-splines with NMI computed from intensities. This is due to the fact that the TPS registration in [START_REF] Mitra | A thin-plate spline based multimodal prostate registration with optimal correspondences[END_REF] is based on control points primarily placed on prostate contour. This results in higher contour accuracy (related to both DSC and HD measures), while B-spline control points are uniformly spread over the prostate image. Since the final aim is accurate localization of biopsy samples and TRE provides a measure of registration accuracy for localized regions, it may be better to choose TRE over DSC or other contour accuracy measures. In this respect our method of B-splines with NMI computed from texture shows statistically significant improvement over B-splines with NMI computed from intensitites; and 1.18 times improvement over TPS registration. However, if the clinical requirement is contour accuracy, TPS registration may be preferred over the B-splines registration.

The proposed method may be extended to 2D-3D registration with the help of an electromagnetic tracking device eliminating the issue of choosing the slice correspondence by an expert. Considering the feasibility of the method of being used for clinical interventions, the entire process may be automated and the algorithm can be implemented on a GPU to be used during real-time interventional needle biopsy procedures.

Fig. 1 A

 1 Fig. 1 A schematic diagram of the proposed algorithm.

Fig. 2

 2 Fig. 2 Even and Odd log-Gabor quadrature filter pairs in spatial domain.

Fig. 3

 3 Fig. 3 Application of quadrature filters on TRUS and MR images. (a) and (c) are the fixed TRUS image and its corresponding quadrature texture; and (b) and (d) are the moving MR image and its corresponding quadrature texture respectively.

7 ±

 7 [START_REF] Kaplan | Real time MRI-ultrasound image guided stereotactic prostate biopsy[END_REF].0 cm3 . The TRUS images were acquired with 6.5MHz side-firing probe with SIEMENS Allegra and TOSHIBA Xario machines and the axial T2 fast relaxation fast spin echo MRI slices with slice thickness of 3 mm, repetition time of 3460 -3860 millisecs and echo time of 113.62 -115.99 millisecs were acquired with GE Signa HDxt. The axial middle slices in TRUS are chosen for which the corresponding axial MR slices are identified by an expert. Ideally the axial MR slice corresponding and parallel to the axial TRUS slice should be obtained by rotating the MR volume in accordance with the rotational angle of the TRUS slice and then resampling the axial MR slice parallel to the axial TRUS slice under observation. However, we have not quantified the TRUS rotational angle in our current experimental process. The prostate is manually segmented from both the moving MR and fixed TRUS images. A NMI-based affine transformation between the TRUS and the MR images is followed by the free-form B-spline deformation. An uniform initial B-spline control grid with a spacing of 64 × 64 is placed on the moving image with an average image size of 256×256. Fig.4shows the uniform B-spline control grids on the moving MR image of size 219×249 at the initial and final resolutions. The Bsplines deform at each resolution to maximize the NMI computed from texture images that are obtained from the magnitudes of quadrature filter responses.

Fig. 4 B

 4 Fig. 4 B-splines control grid with 2 refinements over the initial placement on the moving MR image of size 219 × 249. (a) shows the initial placement of the B-spline grids with 64× 64 pixel spacing on the moving MR image and (b) shows the final set of B-spline control grids on the transformed moving image.

Fig. 5

 5 Fig. 5 Qualitative B-spline registration using NMI from intensity images, NMI from texture images and TPS registration. Patients 6, 8, 10, 11, 13 and 15 in columns. 1 st row shows the fixed TRUS slices, 2 nd row shows the moving MR slices, 3 rd -4 th rows show the fused MR and the checkerboards for B-spline using NMI from intensity images, 5 th -6 th rows show the fused MR and the checkerboards for B-spline using NMI from texture images and 7 th -8 th rows show the fused MR and checkerboard for TPS registration.

Table 1

 1 A comparison of registration accuracies of the B-spline registration with NMI computed from intensities, from texture and TPS registration. µ is the mean and σ is the standard deviation of the measures.

				B-spline Registration				TPS Registration
	Patient#	NMI from Intensities Rueckert et al. [32]	NMI from Texture		Mitra et al. [27]
		DSC	TRE (mm)	TLE (mm)	HD (mm)	DSC	TRE (mm)	TLE (mm)	HD (mm)	DSC	TRE (mm)	TLE (mm)	HD (mm)
	1	0.902 5.07	0.10	9.35	0.896 4.09	0.22	7.31	0.971 9.36	0.22	1.84
	2	0.980 0.37	0.09	1.07	0.964 1.50	0.14	1.30	0.957 3.98	0.10	2.32
	3	0.973 -	-	2.08	0.953 -	-	5.45	0.974 -	-	2.09
	4	0.985 1.91	0.28	1.07	0.962 1.40	0.29	3.96	0.982 5.21	0.49	1.40
	5	0.889 9.08	0.04	8.06	0.869 5.65	0.05	6.76	0.972 2.11	0.07	2.35
	6	0.869 6.11	0.04	6.16	0.975 2.70	0.07	2.23	0.979 1.17	0.05	2.32
	7	0.959 0.90	0.12	5.04	0.889 1.94	0.09	8.82	0.977 4.43	0.12	2.61
	8	0.976 4.70	0.03	3.38	0.964 1.11	0.10	5.81	0.978 3.57	0.05	2.96
	9	0.960 -	-	4.93	0.982 -	-	1.32	0.978 -	-	2.39
	10	0.952 8.29	0.09	5.98	0.981 4.04	0.13	1.04	0.972 6.09	0.04	1.98
	11	0.962 6.12	0.04	3.03	0.950 1.69	0.04	3.72	0.972 2.98	0.12	2.22
	12	0.944 1.58	0.25	4.60	0.934 0.31	0.16	5.18	0.971 2.44	0.12	5.00
	13	0.961 1.00	0.05	4.83	0.878 3.42	0.05	7.39	0.980 3.06	0.07	1.84
	14	0.896 7.32	0.05	6.50	0.965 1.85	0.05	3.12	0.986 1.75	0.07	0.82
	15	0.942 5.01	0.05	6.06	0.948 2.32	0.13	5.02	0.968 2.29	0.07	2.22
	16	0.974 7.27	0.13	2.90	0.950 4.71	0.13	6.50	0.970 1.86	0.07	2.71
	17	0.894 4.32	0.12	7.18	0.872 3.22	0.23	7.64	0.982 0.18	0.32	1.04
	18	0.985 0.51	0.26	2.32	0.975 1.77	0.29	2.87	0.982 0.91	0.26	1.64
	19	0.936 5.12	0.19	7.01	0.969 2.91	0.43	7.09	0.983 1.47	0.23	1.66
	20	0.939 5.03	0.20	5.51	0.975 2.83	0.21	2.45	0.973 3.11	0.30	2.08
	µ	0.944 4.43	0.12	4.85	0.943 2.64	0.16	4.75	0.975 3.11	0.15	2.17
	σ	0.036 2.77	0.08	2.30	0.039 1.37	0.11	2.40	0.007 2.18	0.12	0.85

ACKNOWLEDGMENTS

This research is a part of the PROSCAN project of the VICOROB laboratory of University of Girona, Catalunya, Spain. The authors would like to thank VALTEC 08-1-0039 of Generalitat de Catalunya, Spain and Conseil Régional de Bourgogne, France for funding this research.


Prof F Javier

Luque

Ph.D.; Josep Rayo

Ph.D.; Lourdes Muñoz

Gloria Rosell

Bruce D Hammock

Ph.D.; Angel Guerrero

Javier F Luque

Ph.D Ramon Pouplana
email: rpouplana@ub.edu

F Javier Luque
email: fjluque@ub.edu

Reactivity versus Steric Effects in Fluorinated Ketones as Esterase Inhibitors: A Quantum Mechanical and Molecular Dynamics Study

Keywords: Fluorinated ketones, Esterase, Quantum Mechanical Computations, Molecular Dynamics, Structure-Based Drug Fluorinated ketones  Esterase  Quantum Mechanical Computations  Molecular Dynamics  Structure-Based Drug Design

Carboxylesterases (CEs) are a family of ubiquitous enzymes with broad substrate specificity, and their inhibition may have important implications in pharmaceutical and agrochemical fields. One of the most potent inhibitors both for mammalian and insect CEs are trifluoromethyl ketones (TFMKs), but the mechanism of action of these chemicals is not completely understood. This study examines the balance between reactivity versus steric effects in modulating the activity against human carboxylesterase 1. The intrinsic reactivity of the ketone moiety is determined from quantum mechanical computations, which combine gas phase B3LYP calculations with hydration free energies estimated with the IEF/MST model. In addition, docking and molecular dynamics simulations are used to explore the binding mode of the inhibitors along the deep gorge that delineates the binding site. The results point out that the activity largely depends on the nature of the fluorinated ketone, since the activity is modulated by the balance between the intrinsic electrophilicity of the carbonyl carbon atom and the ratio between keto and hydrate forms. However, the results also suggest that the correct alignment of the alkyl chain in the binding site can exert a large influence on the inhibitory activity, as this effect seems to override the intrinsic reactivity features of the fluorinated ketone. Overall, the results sustain a subtle balance between reactivity and steric effects in modulating the inhibitory activity of TFMK inhibitors.

Following the reviewer's suggestion, we have made a thorough revision in order to clarify her/his points (see below), to make a more comprehensive analysis of the results and to clarify the main findings in Conclusions.

a) The sources of chemicals has been included in Methods providing appropriate citations and describing the synthetic protocol in one case (compound 4). Chemical characterization of this compound is given in Supporting Information.

b) The text provides more details about the protocol used in docking calculations. Starting structures of the TFMK inhibitors were generated in an extended geometry, which was subsequently refined through energy minimization with MMF94 force field. The adoption of an extended geometry as a starting guess seems reasonable, as the alkyl chain of palmitic acid is rather elongated in the complex with hCE1 (the distance from the carbon bearing the carboxyl group to the last methylene in the conformation found in 2DQY is 11.6 Å, which compares with a value of 16.4 Å for the same distance in a fully extended conformation). The conformational flexibility of the ligand was explored in docking calculations. However, two constraints were imposed to mimic the covalently-bound tetrahedral intermediate resulting from the nucleophilic attack of the hydroxyl group to Ser221 to the carbonyl carbon: i) the covalent bond between Ser221 (oxygen) and the carbonyl group (C), and ii) the hydrogen bonding betwen the carbonyl oxygen of the inhibitor and the backbone NH units of Gly142 and Gly143. These contraints also permit to enhance the conformational sampling in exploring the alignment of the alkyl chain. c) A more detailed description of the parametrization of the tetrahedral intermediate formed between the TFMK inhibitor and the catalytic Ser221 is now provided. Particular attention was paid to the charges of the (Ser)-O-C(O)-(TFMK) unit, as this is relevant for the proper description of the hydrogen bond interactions with the NH units of Gly142 and Gly143.

d) The description of the IEF/MST formalism has been shortened. e) We have introduced a more detailed discussion of present results and their interpretation in light of the data given in ref. [START_REF] Wheelock | [END_REF] (26 in the original version). We agree with the reviewer that the experimental work reported in ref. [START_REF] Wheelock | [END_REF] strongly supports the hypothesis that the inhibitory potency is related to the degree of ketone hydration. However, this is not reflected in the regression equations between inhibitory potency and energetic properties, since they all rely on linear equations. Our results, in contrast, demonstrate a parabolic behavior, which can then be interpreted from the balance between the degree of ketone hydration and the intrinsic ketone reactivity. The point we want to stress is that our study reconciles theoretical predictions and experimental data and provides a clear support to the interpretation about the role played by ketone hydration on the inhibitory activity of TFMKs. f) Figures 6, 7 and 9 have been redrawn following the reviewer's suggestions.

g) The abbreviation SE is explained in Tables 2,3.

Introduction

Carboxylesterases (CEs; EC 3.1.1.1) include a heterogeneous group of isozymes that hydrolyze a wide range of aliphatic and aromatic esters, amides and thioesters [1,[START_REF] Quinn | (ed) Enzymes, Enzyme Mechanisms, Proteins and Aspects of NO Chemistry[END_REF]. They comprise a multigene family that covers key hydrolytic enzymes with extremely broad substrate specificity [START_REF] Satoh | [END_REF]4]. In addition, CEs are important in pest resistance since they participate in detoxification processes of carbamate and organophosphate insecticides, and insect strains resistant to these compounds have shown high levels of CEs [5,6]. Therefore, the development of new CE inhibitors may lead to drugs with therapeutic and/or agrochemical interest [7,8]. Esterases are important also in the degradation of a large variety of insect pheromones and hormones, and thus are involved in many chemical and biochemical processes.

An important group of chemicals that promote CE inhibition are fluorinated ketones, which inhibit a variety of esterases such as acetylcholinesterase [9], juvenile hormone esterase (JHE) [10] and human liver microsomal CEs [11,12], as well as enzymes such as chymotrypsin and trypsin [13,14] and enzymes that metabolize chemical mediators including fatty acid amide hydrolase and diacyl glycerol [15,16]. Trifluoromethyl ketones (TFMKs) also reversibly inhibit pheromone-degrading esterases in male olfactory tissues [17][18][19]. Therefore, these compounds represent a potential approach for the disruption of pheromone communication as a new strategy for pest control [20,21]. Nevertheless, further research is required to increase the clinical and agricultural application of these compounds [22].

TFMKs behave as transition-state analogues, since the enzyme inhibition involves formation of a tetrahedral adduct between the serine residue present at the active site and the highly electrophilic carbonyl moiety [23,24]. Support for this mechanism of action comes from Xray crystallographic structures of the covalent complex of porcine pancreatic elastase and a peptidyl ,-difluoro--ketoamide [25], and of the JHE from the tobacco hornworm Manduca sexta with 3-octyl-1,1,1-trifluoropropan-2-one [26]. Structure-activity relationship studies indicate that the potency of TFMK inhibitors is modulated by lipophilicity, with the optimal activity being associated with intermediate lipophilicity values (3<log P<5), molar refractivity and the degree of fluorination [12,[START_REF] Székács | Rational Approaches to Structure, Activity and Ecotoxicology of Agrochemicals[END_REF]. These findings agree with the high hydrophobicity of the gorge that leads to the catalytic site in CEs, and suggest a direct influence of substituents attached to the ketone moiety on the chemical reactivity towards formation of the tetrahedral adduct with the enzyme. The inhibitory potency has also been related to the degree of ketone hydration [START_REF] Wheelock | [END_REF], so that inhibitors that favour the tetrahedral geometry of the hydrate (gem-diol) should stabilize the enzyme-inhibitor complex, thereby resulting in more potent inhibitors. Finally, it has also been suggested that in -thiosubstituted TFMKs one of the hydroxyl groups in the gem-diol might form an intramolecular hydrogen bond with the sulfur atom, and that the strength of such hydrogen bond shows a good correlation with the inhibitory potency of -substituted TFMKs involving several groups (O, S, SO and SO 2 ) [29].

In order to gain further insight into the precise nature of the molecular determinants that mediate the CE inhibitory activity of TFMKs, this work examines the contribution of both reactivity and steric properties of a series of TFMK inhibitors. The activity of the compounds has been determined experimentally against human carboxylesterase 1 (hCE1). The reactivity properties of the fluorinated ketone moieties are examined by means of density functional computations in conjunction with continuum solvation computations, paying particular attention to the balance between keto and hydrate forms of the ketone moiety. Moreover, the effect of the alkyl chain in TFMKs on the inhibitory potency is investigated by combining docking calculations with molecular simulations. The results are then used to discuss the relationship between the molecular structure of TFMKs and the hCE1 inhibitory potency.

Methods

Chemicals

Compounds 1 and 5 were prepared as already reported by us [30,31], and the synthesis of compounds 2 [12] and 3 [32] has been already described. Compounds 6 and 7 were obtained from n-decanal and n-octanamide, respectively (manuscript in preparation). Difluoromethyl ketone 8 was prepared as previously described [33] and 2-dodecanone (9) was commercially available from Acros Organics (Geel, Belgium). Compound 10 was prepared by reaction of the corresponding thiol with 4-chloro-1,1,1-trifluoro-2-butanone (manuscript in preparation).

Compounds 11-13 were taken from ref [34]. TFMKs 14, 16 and 17 were prepared by transmetallation reaction of the corresponding iododerivatives with tert-BuLi followed by treatment with ethyl trifluoroacetate [31]. Compound 15 was synthesized from 1-octadecene with trifluoroacetaldehyde ethyl hemiacetal in the presence of a Lewis acid followed by Dess-Martin oxidation [35] of the resulting trifluoromethyl carbinol (submitted for publication).

Finally, methylketone 18 was obtained by methyllithium addition to the corresponding carboxylic acid [36].

Synthesis of 1,1,1-trifluoro-3-octyloxypropan-2-one (4) BF 3 •Et 2 O (5.5 L, 0.044 mmol) was slowly added to a stirred solution of 1-octanol (400 mg, 3.07 mmol) and 1,1,1-trifluoroepoxypropane (280 mg, 2.5 mmol) in CH 2 Cl 2 (1 mL). The reaction mixture was stirred at 40-45 ºC for 16 h. The reaction mixture was partitioned between CH 2 Cl 2 (10 mL) and H 2 O (10 mL). The organic solution was washed with H 2 O (2 x 10 mL), dried (MgSO 4 ), filtered and evaporated to afford a yellow oil, which was purified by column chromatography on silica gel eluting with hexane:ether 9:1. The corresponding trifluoromethyl carbinol (386 mg, 63%) was thus obtained as colorless oil. A 15% w/w Dess-Martin reagent in CH 2 Cl 2 (24 mL, 3.68 g, 1.05 mmol) was then added to a solution of the carbinol (188 mg, 0.77 mmol) in CH 2 Cl 2 (0.2 mL). The reaction mixture was stirred at room temperature for 19 h. The resultant cloudy solution was diluted with Et 2 O (10 mL), washed with 10% w/v aq. Na 2 S 2 O 3 (10 mL), NaHCO 3 saturated solution, and brine, and dried over MgSO 4 . After filtration and concentration, the residual oil was purified by flash column chromatography eluting with hexane-ether 85:15 to afford compound 4 (146 mg, 78%) as colorless oil.

Biological data

Assays were performed in 96-well polystyrene microtiter flat bottom microplates (Maxisorp, Labclinics, Barcelona, Spain) at 32ºC using a 0.4 mM solution of p-nitrophenyl acetate (PNPA) in acetonitrile as substrate. Stock solutions of the possible inhibitors were prepared in ethanol from which 2 µL were used. For the assays, 50 µL of the enzyme in sodium phosphate buffer (containing 10 µg of hCE1) with the inhibitor was incubated for 15 min, and then 10 µL of the substrate solution was added for a final concentration of 0.4 mM. Each solution was diluted with sodium phosphate buffer (pH 7.8, 100 mM) to a final volume of 200 µL. The rate of PNPA hydrolysis was monitored every 30 s for 5 min in a Spectramax Plus plate reader (Molecular Devices, Sunnyvale, CA) at 405 nm (p-nitrophenyl alcohol absorption). Each inhibitor concentration was tested in triplicate, and the assay was repeated 3 times along with positive (in the presence of enzyme but no inhibitor) and negative (no enzyme, no inhibitor) controls. Inhibitor concentrations were adjusted so that the linear region of the curve had at least four points above and below the IC 50 . The IC 50 values were determined by regression analysis of at least eight points in the linear region of the curve. For the specific case of compounds 14-18, the experimental procedure described by Nishi et al. was adopted [37].

Quantum mechanical computations

Theoretical calculations were performed for model compounds built up by replacing the alkyl chain of TFMKs by an ethyl group (Fig. 1) in order to explore the reactivity properties of the carbonyl moiety present in TFMKs. This choice is justified not only by the reduction in cost of the quantum mechanical computations, but also to establish the relationship between the reactive properties of the substituted ketone moiety and the inhibitory activity. The geometries of model compounds were fully optimized at the B3LYP/6-31G(d) level [38,39].

The nature of the stationary points was verified by inspection of the vibrational frequencies, which were used to calculate zero point, thermal and entropic corrections within the framework of the harmonic oscillator-rigid rotor at 1 atm and 298 K. These corrections were added to the electronic energies to estimate the free energy differences (

gas G 
) in the gas phase.

The relative stabilities in solution were estimated by combining the free energy difference in the gas phase and the differences in solvation free energy (

sol G 
) determined by using the IEF/MST solvation model [40,41], which relies on the Integral Equation formalism (IEF) of the polarizable continuum method [42,43]. The IEF-MST model computes

sol G  from the addition of electrostatic ( ele G  ), cavitation ( cav G 
) and van der Waals (

vW G  ) components. ele G 
is determined from the interaction between the charge distribution of the solute and the electrostatic response of the solvent, which is treated by a set of point charges spread over the surface of the cavity that separates solute and solvent.

cav G 
is computed following Claverie-Pierotti's scaled particle theory [44,45]. Finally,

vW G 
is computed using a linear relationship to the solvent-exposed surface of each atom [41,46]. IEF/MST computations were performed using the B3LYP/6-31G(d) optimized version of the MST(IEF) model [41]. Gas phase calculations were carried out using Gaussian-03 [START_REF] Frisch | Gaussian03[END_REF], and IEF/MST calculations were performed using a locally modified version of this program.

Molecular modeling

Docking calculations were used in conjunction with molecular dynamics simulations in order to examine the alignment of the alkyl chain along the gorge and its effect on the activity of selected TFMK inhibitors. To this end, the X-ray crystallographic structures of hCE1 complexed with palmitic acid (PDB entry 2DQY; solved at 3.0 Å resolution [48]) and benzoic acid (PDB entry 1YAJ; solved at 3.2 Å resolution [49]) were used in the docking study. The former was chosen due to the fact that the alkyl chain of palmitic acid delineates the gorge leading to the active site, whereas the latter was selected because the benzoic acid is found forming a covalently-bound tetrahedral intermediate with the catalytic serine Ser221 in subunits C, F and J (in the rest of subunits the benzoic acid is unbound in the catalytic site).

Accordingly, these structures provide valuable information for the proper positioning of the covalent adduct formed by the TFMK derivatives investigated here with the hCE1 enzyme.

Docking studies were carried out using Gold 4.1 (CCDC, Cambridge) and the Goldscore scoring function [50][51][52]. The X-ray structure of hCE1 was used as a template for docking calculations after removal of ligands, ions and water molecules. The definition of the binding site to be considered for docking of TFMK inhibitors was made taking advantage of the known spatial arrangement of both palmitic acid and benzoic acid in the X-ray structures 2DQY and 1YAJ. The structure of each inhibitor was initially built up with the alkyl chain in an extended conformation using MOE (Chemical Computing Group, Montreal), and the geometry was subsequently refined by energy minimization using the MMFF94s [53] force field. This extended conformation seemed adequate as starting geometry, as the alkyl chain of palmitic acid is rather elongated in the complex with hCE1 (the distance from the carbon bearing the carboxyl group to the last methylene in the conformation found in 2DQY is 11.6 Å, which compares with a value of 16.4 Å for the same distance in a fully extended conformation). It is worth noting, however, that whereas the protein was kept rigid, Gold accounts for the conformational flexibility of the ligand around rotatable bonds during docking calculations. In order to explore the alignment of the inhibitor in a suitable orientation that mimics the covalently-bound tetrahedral intermediate arising from the nucleophilic attack of the hydroxyl group of Ser221 to the carbonyl unit of the TFMK inhibitor, a covalent bond was imposed between the oxygen atom of Ser221 and the carbonyl carbon atom in the inhibitor. In addition, a hydrogen bonding constraint was imposed to bias the orientation of the carbonyl oxygen in the tetrahedral intermediate towards the backbone NH units of Gly142 and Gly143. Thus, this constraint should account for the stabilization played by those residues upon formation of the tetrahedral adduct in the enzyme-substrate complex. Overall, the application of these constraints permits to enhance the conformational search in sampling the alignment of the alkyl chain in the binding site.

Refinement of the docked complexes was accomplished by means of molecular dynamics simulations using the parmm99 force field of the Amber-9 package [START_REF] Case | AMBER, Version 9. University of California[END_REF]. The general Amber force field (gaff; [55,56]) was used to assign parameters for the covalently-bound adduct formed between the ketone moiety of the TFMK inhibitor and the hydroxyl group of Ser221.

In order to refine the gaff parameters, the adduct was modeled by replacing the alkyl chain by an ethyl group and by adding acetyl and methyl groups to the N-and C-terminals of Ser221.

The geometry of this model compound was subsequently optimized at the HF/6-31G(d) level, though the torsional angles of the backbone were kept frozen in order to avoid an unrealistic rearrangement in the molecular geometry of the peptide units. The optimized geometrical parameters were then used to refine the bond lengths and angles of the tetrahedral intermediate. Moreover, in order to improve the description of the charge distribution of the (Ser)-O-C(O)-(TFMK) unit and to describe more accurately the hydrogen-bonding stabilization arising from the interaction with the NH units of Gly142 and Gly143, the partial atomic charges were slightly adjusted from the RESP charges [57] determined at the HF-6-31G(d) level (note that this level of theory was already adopted in the derivation of the gaff force field).

The standard ionization state at neutral pH was taken for all the ionizable residues of the enzyme but His468, which was protonated. The complexes were then immersed in a preequilibrated box (ca. 109 3 Å 3 ) containing the inhibitor-enzyme complex and around 17000 TIP3P [58] water molecules. Six Na + ions were added using the CMIP strategy to neutralize the system [59], which was then energy minimized and equilibrated using a multistep protocol. First, the location of hydrogen atoms was energy minimized for 2000 steps of steepest descent. Next, the positions of water molecules and ions were relaxed for 2000 steps of steepest descent plus 5000 steps of conjugate gradient. Finally, the whole simulated system was refined via energy minimization by combining 2000 steps of steepest descent plus 8000 steps of conjugate gradient. At this point, the system was thermalized at 298 K by running four 25 ps molecular dynamics simulations where the temperature was increased from 100 to 150, 200, 250 and 298 K, respectively. Subsequently, a 10 ns MD simulation was carried out.

The system was simulated in the NPT ensemble using periodic boundary conditions and Ewald sums for treating long-range electrostatic interactions (with the default Amber-9 parameters). The structural analysis was performed using in house software and standard codes (PTRAJ module) of Amber-9.

Results and discussion

TFMK inhibitors

Two series of compounds were considered to investigate the relationship between the structure and the inhibitory potency of TFMKs. The first series includes thirteen inhibitors that combine a short alkyl chain with a ketone moiety, which is the main source of chemical diversity due to factors such as the presence or absence of fluorine atoms in the terminal methyl group and chemical modifications in positions vicinal to the carbonyl group (compounds 1-13; Fig. 1). In contrast, the second series includes inhibitors that share the same trifluoromethyl ketone moiety, but differ in the nature of the long alkyl chain (compounds 14-17; Fig. 2). In the former case, it is reasonable to expect that the high flexibility of the short chain should not impede the proper alignment of the inhibitor in the catalytic site. Accordingly, the differences in inhibitory activity should mainly arise from the intrinsic reactivity of the ketone moiety, which in turn should influence the formation of the tetrahedral adduct with the catalytic serine. In the second case, the inhibitory potencies are expected to reflect differences in the proper arrangement of the long alkyl chain, particularly due to the presence of the double bond in positions  and  (relative to the carbonyl group) in compounds 15 and 16, respectively, compared to the activity of 14.

Reactivity studies

The efficiency of TFMKs as inhibitors of CEs is related to the susceptibility of the ketone moiety to participate in the nucleophilic attack by the catalytic serine in the enzyme active site, thus leading to a tetrahedral adduct that mimics the transition state of the enzymesubstrate complex. Accordingly, in the absence of other factors that might prevent a correct positioning of the ketone moiety in the catalytic site, those inhibitors that favour the formation of the tetrahedral intermediate should a priori exhibit the largest biological effect. Under this rationale, we have examined the equilibrium between keto and hydrate species of model compounds 1-13 to gain insight into the reactivity properties of the ketone moiety. To this end, we have adopted the thermodynamic cycle shown in Fig. 3, which allows us to combine the intrinsic (gas phase) stability of the keto and hydrate species with the hydration preferences of the two species. For those compounds able to form an intramolecular hydrogen bond, a limited conformational search was also performed both in the gas phase and in aqueous solution in order to examine the relative stability of conformers with and without intramolecular hydrogen bonds in both phases.

The stability of the hydrate form in the gas phase ( gas G 

; Table 1) depends on the number of fluorine atoms positioned  to the carbonyl group. Thus, the structure with the lowest stability ; Table 1) upon solvation in aqueous solution. This effect can be attributed to the weakening of the intramolecular hydrogen bond upon solvation in water. In contrast, the relative stability of the keto and hydrate forms is less affected for those compounds whose molecular structure does not permit the formation of such intramolecular hydrogen bond (compounds 5, 8, 9 and 10).

The preceding discussion points out that the balance between the hydrate and keto species depends on a subtle balance of different factors depending on the chemical nature of the groups introduced in the ketone moiety. The degree of fluorination of the methyl group vicinal to the carbonyl unit influences the electron density on the carbonyl carbon, which becomes more susceptible to the nucleophilic attack of a water molecule as the number of fluorine atoms increases. In the gas phase the stability of the hydrate form is also enhanced by the formation of an intramolecular hydrogen bond with a neighbouring group suitably placed in the ketone moiety to act as hydrogen-bond acceptor. The magnitude of this latter effect can however be largely counterbalanced upon solvation in polar solvents due to the competing interactions with polar solvent molecules. The combined effect of these factors is reflected in a sizable difference in the preference between keto and hydrate species for the series of compounds considered in this study, as can be seen in Table 1.

It is well known that the hydration constant of carbonyl compounds in aqueous solution can vary in a wide range (10 -6 to 10 3 ) [60]. The accurate determination of the ratio between keto and hydrate forms is difficult since the rate of hydration is usually too fast to allow determination of the equilibrium constant by isolating one of the forms. Moreover, in many cases the amount of keto or hydrate forms may be very small, only detectable by very sensitive spectroscopic methods. On the basis of these considerations, we attempted to measure experimentally the ratio between keto and hydrate forms for few selected compounds using 19 F NMR in chloroform solution (note that these assays were performed for compounds with R= hexyl). Under these conditions, the results indicate that hydration of the keto form follows the order sulfone (3; 100% hydrate)  sulfoxide (2; 90-96% hydrate) > ether (4; 78-80% hydrate) > thioether (1; 30-50% hydrate) > unsubstituted TFMK (5; 0% hydrate)  difluoromethyl ketone (8; 0% hydrate). Keeping in mind the different polarity of the solvent used in computations (Table 1) and in experimental measurements, one can conclude that there is satisfactory qualitative agreement between the theoretically predicted and experimentally determined preferences between the keto and hydrate forms of the studied compounds.

In order to rationalize the hydration extent of compounds 1-13, the free energy difference between keto and hydrate species was compared with the partial charge of the carbonyl carbon, which was determined by using the atoms-in-molecules theory [START_REF] Bader | Atoms in Molecules: A Quantum Theory[END_REF]. Fig. 4 shows that there is a correlation between the free energy difference determined in aqueous solution and the charge of the carbonyl carbon (the only exception to this behaviour is compound 6, which likely reflects the influence of the vicinal carbonyl group). Clearly, this trend can be realized from the influence exerted by the degree of fluorination in  position and the substitution by electron-withdrawing groups at  position of the carbonyl group on the carbon charge, which becomes more electrophilic and therefore more prone to hydration.

The inhibitory potency determined experimentally for the TFMK inhibitors that incorporate the ketone moieties 1-13 is given in Table 2. In general, there is a close similarity between the experimental activities for the two series of compounds reported in Table 2, which differ in the length of the short alkyl chain (R is hexyl in this work, and butyl in ref. [34]). Such an agreement supports the assumption that the differences in inhibitory potency are mainly due to the intrinsic reactivity properties of the ketone moiety of TFMKs.

Previous studies have examined the dependence on the inhibitory potency of TFMKs and the relative stability between ketone and hydrate species [START_REF] Wheelock | [END_REF]. For a small series of six compounds, which includes sulfide (1,11), sulfone [START_REF] Satoh | [END_REF]12) and ether (4, 13) methyl ketone moieties and their trifluorinated derivatives, a linear dependence was found between the inhibitory activity against human carboxylesterase (from human liver microsomes) and the gas phase electronic energy difference between ketone and hydrate species [START_REF] Wheelock | [END_REF]. Similar linear dependencies were also found with the gas phase free energy difference and with the relative stability estimated in aqueous solution using two solvation models [START_REF] Wheelock | [END_REF]. These findings mainly reflect the activating effect due to the fluorination of the terminal methyl group in the ketone. Nevertheless, the linear dependence between the inhibitory activity and the energetic difference of keto and hydrate species is counterintuitive, as it would a priori suggest an enhancement in the inhibitory potency as the stability of the hydrate form is enlarged. In fact, the kinetic study reported for the trifluorinated sulfoxide derivative 2, whose predicted activity (pIC 50 = 10.00) was much larger than the experimental value (pIC 50 = 6.51), suggested that the dehydration of the gem-diol form was likely the rate-limiting step of binding to the enzyme.

In contrast to the preceding findings, comparison of the inhibitory potencies determined experimentally for compounds 1-13 with the free energy difference reported here for keto and hydrate species in aqueous solution reveals the existence of a parabolic dependence of the IC 50 values (Fig. 5). Compared to the theoretical linear relationships reported in [START_REF] Wheelock | [END_REF], the distinct functional dependence found in this study can be attributed to the inclusion of a larger, more diverse set of compounds. The parabolic relationship between inhibitory potencies and keto-hydrate relative stabilities can be understood from the balance between two opposite factors: the hydration of the keto species and its intrinsic reactivity against a serine residue of the enzyme. Thus, an increase in the electrophilicity of the carbonyl carbon should enhance the reactivity for the attack of the catalytic serine residue, but at the same time it also implies a displacement of the compound towards the hydrate form in aqueous solution.

Accordingly, the activity of those ketone moieties poorly activated for a nucleophilic attack will be limited by the intrinsic reactivity, though they will tend to populate the keto species in solution. Likewise, those compounds containing a ketone moiety with enhanced reactivity towards a nucleophile will be moderate or even poor inhibitors because the concentration of the inhibitor will mainly populate the hydrate species, and the inhibitory activity will therefore be largely determined by the dehydration of the gem-diol form to generate the ketone species.

Overall, Fig. 5 allows us to reconcile the predictions made from theoretical calculations for the balance between keto and hydrate species of TFMK compounds with the experimental evidences reported in [START_REF] Wheelock | [END_REF], which support the crucial role played by the degree of ketone hydration on the inhibitory potencies. Moreover, it can also be deduced that if the ketohydrate conversion is sufficiently dynamic as to ensure the equilibrium between those species, increased incubation times would then result in enhanced inhibitory potencies, because dehydration of the gem-diol species would yield the ketone species, as found in previous experimental studies [12].

Modeling studies

The preceding results reflect the influence of the electronic properties exerted by substituents on the reactivity of the ketone moiety. In particular, the results point out the important role due to fluorination in enhancing the reactivity of the compounds relative to their nonfluorinated partners. This trend would suffice to justify the drastic reduction in the hCE1 inhibitory activity measured for compound 18 relative to TFMKs 14 and 17 (see Fig. 2 and Table 3). Such a reduction cannot be attributed to the different length of the alkyl chains Therefore, it can be hypothesized that the conformational restraints imposed by the double bond could affect the proper alignment of the TFMK inhibitor in the catalytic site, thus preventing the formation of the tetrahedral intermediate with the catalytic serine residue.

To corroborate this hypothesis, docking calculations were combined with molecular dynamics simulations in order to examine the ability of the catalytic site to accommodate the covalently bound adducts generated from inhibitors 14-16. Inspection of the docked poses revealed that the alkyl chain was generally suitably aligned along the hydrophobic gorge. However, due to the pharmacophore constraints imposed to the carbonyl group of the ketone in docking computations (see Methods), two main orientations were found for the terminal trifluorinated methyl group, as noted in Fig. 7 for compound 14 (similar results were obtained for 15 and 16; data not shown). In the two poses the orientation of the trifluoromethyl ketone moiety permits the correct positioning of the hydroxyl oxygen of Ser221 for the nucleophilic attack to the carbonyl carbon. Moreover, the hydroxyl unit in Ser221 is assisted by a hydrogen bond with the imidazole ring of His468. Nevertheless, the two poses differ in the orientation of the trifluoromethyl unit in the binding site. In one case (shown in orange in Fig. 7) no steric clash is observed between the trifluoromethyl group and the neighbouring residues in the catalytic site. However, in the other case (shown in yellow in Fig. 7) the trifluorinated methyl is positioned close to the NH unit of Gly142 and the imidazole ring of His468, suggesting the existence of unfavorable clashes. Accordingly, this latter pose was not considered for further analysis in molecular dynamics simulations.

The structural integrity of the most favoured pose for the tetrahedral adduct formed between inhibitors 14-16 and the catalytic serine of hCE1 was examined by means of 10 ns molecular dynamics simulations. Inspection of the time evolution of the potential energy of the simulated systems supported the stability of the trajectories in the last 5 ns (see Fig. S1 in Supporting Information). The profile for the positional root-mean square deviation (RMSD)

of the backbone atoms ranged from 1.8 to 2.0 Å and remained stable during the last 3 ns, thus reflecting a small relaxation of the overall protein skeleton compared to the X-ray crystallographic structure (see Fig. S1). The RMSD of the heavy atoms for the residues that delineate the walls of the catalytic site and the gorge was smaller (ranging from 1.3 to 1.8 Å)

and remained stable along the last 4 ns of the trajectories (Fig. 8).

Though the preceding results point out that the overall stability of the trajectories was similar for the three hCE1-inhibitor complexes, careful inspection of the residues at the binding site reveals a significant difference in the conformation of His468. This residue plays a relevant functional role, as it contributes to enhance the nucleophilicity of Ser221 through the formation of a hydrogen bond, and in turn interacts with Glu354. Thus, inspection of the three subunits present in the X-ray crystallographic structure 2DQY indicates that the distances from O(Ser221) to N  (His468) vary between 2.9 and 3.2 Å, and the distances from N  (His468) to O(Glu354) range from 2.6 to 2.8 Å. Interestingly, the triad formed by Ser221, His468 and Glu354 is structurally stable along the trajectories run for the adducts generated from TFMKs 14 and 16 (see Fig. 9). Thus, the hydrogen bond between Ser221 and His468 remains stable along the whole trajectory (with an average distance close to 2.9 Å; see Fig. 10). In contrast, the structural arrangement of the catalytic triad is completely lost in the simulation run for inhibitor 15 (Fig. 9), as noted in the large interatomic distance between O(Ser221) and N(His468) (> 6 Å; Fig. 10).

The structural alteration of the catalytic triad, which is found since the beginning of the production run, would imply the loss of the catalytic power.

Final remarks

Previous studies have indicated that an increased hydrophobicity correlates with inhibitor potency, so that compounds containing longer, more hydrophobic alkyl chains are more potent inhibitors of CEs [12,[START_REF] Székács | Rational Approaches to Structure, Activity and Ecotoxicology of Agrochemicals[END_REF]. Clearly, this trend reflects the lipophilic nature of the gorge, as found in the X-ray crystallographic structure of hCE1, where the binding site is lined by residues such as Ala93, Leu97, Val146, Val254, Leu255, Leu299, Phe303, Leu304, Ile359, Leu363, and Met425. However, other factors are expected to play a decisive role on the inhibitory activity of TFMKs.

The results reported here allow us to stress the relevant role played by substituents present in the ketone moiety, which modulate the intrinsic reactivity against a nucleophilic attack, but also the balance between keto and hydrate species. Thus, the parabolic dependence found between the inhibitory potency and the relative stability between keto and hydrate forms of TFMKs strongly support the hypothesis that the inhibitory potency is related to the degree of ketone hydration, as noted in previous studies [START_REF] Wheelock | [END_REF]. Such dependence suggests a direct implication of the electrophilicity of the carbonyl carbon atom of the ketone moiety in two counterbalancing effects: the susceptibility for the nucleophilic attack of the catalytic serine residue of the enzyme, and the degree of hydration of the ketone species.

In this context, present findings support the assumption of the keto species as the bioactive form of TFMK inhibitors [START_REF] Wheelock | [END_REF] provided that the TFMK could fit properly in the catalytic cavity. The presence of apparently minor chemical modifications in positions close to the ketone moiety may thus have an unexpected effect on the inhibitory potency. This is the case of the positional isomerism associated with the displacement of the double bond from  to  position in compounds 15 and 16, which gives rise to a drastic change in the inhibitory activity. In this case, the intrinsic reactivity properties of the trifluorinated methyl ketone moiety are overridden by the steric factors associated with accommodation of the inhibitor. In turn, these findings pave the way to take advantage of structural differences in the binding site in order to modulate the selectivity of this class of compounds towards different CEs.

Understanding the molecular determinants that modulate the interaction of TFMK inhibitors with CEs is necessary in order to develop compounds with improved pharmacological and agricultural properties. Based on the preceding findings and their implication for the mechanism of action of TFMKs, we feel that selectivity of inhibitors with a central pharmacophore of a polarized ketone for different target enzymes can be modulated through the inclusion of suitable substituents in the chemical structure that contains the ketone moiety.

Finally, it is likely that general conclusions drawn from the interaction of TFMK inhibitors and CEs can be extended to the design of inhibitors of other targets such as fatty acid amide hydrolase [15] and diacyl glycerol [16].

Tables

Figure captions

5 . 1 (

 51 of the hydrate species corresponds to the non-fluorinated ketone 9 and introduction of three fluorine atoms in the terminal methyl group (5) stabilizes the hydrate form by 8.3 kcal mol -1 . This effect is also observed in the comparison of the relative stabilities for pairs of compounds 1 and 11, 3 and 12, and 4 and 13, where fluorination largely stabilizes the hydrate form compared to the keto species. Similarly, formation of an intramolecular hydrogen bond between one of the hydroxyl groups of the hydrate and the heteroatom in  position to the carbonyl group also enhances the stability of the hydrate, as noted in the comparison of thioether (1), sulfoxide (2), sulfone (3) or ether (4) derivatives with the unmodified compound The stability of the hydrate species is, nevertheless, reduced by around 2-5 kcal mol -

(C 18 H

 18 35 in 14, and C 16 H 31 in 17 and 18), since the terminal part of the long alkyl chain should be easily accommodated in the mouth of the gorge, as noted upon inspection of the different orientations found for the terminal part of the chain of palmitic acid (C 15 H 31 COOH) in the three subunits of the X-ray crystallographic structure (PDB entry 2DQY; Fig.6). Rather, the enhanced inhibitory potency of compounds 14 (IC 50 = 14 nM) and 17 (IC 50 = 9.9 nM) compared to 18 (IC 50 > 1000 nM) can be explained by the enhanced reactivity of TFMK associated with the presence of fluorine atoms in  position to the carbonyl. Surprisingly, the data in Table3also shows a drastic reduction in the inhibitory potency of compound 15 (IC 50 > 1000 nM), which compares with the similar change in activity observed for the non-fluorinated compound 18. This unexpected trend cannot be realized from electronic factors, as the trifluorinated methyl ketone moiety is preserved in compound 15, which should therefore have similar propensities for the nucleophilic attack by the catalytic Ser221. Likewise, the weak inhibitory activity of 15 can hardly be attributed to the hydrophobicity of the compound, as the length of the alkyl chain is identical in compounds 14-16, which only differ in the number and location of the double bonds present in the chain.Inspection of the chemical structures of those compounds suggests that the poor inhibitory potency of 15 could be attributed to the precise location of the double bond between  and  positions in the alkyl chain, as the inhibitory activity of compound 16 (IC 50 = 11.8 nM), where such a double bond involves carbon atoms  and , nearly matches the activity of 14 and 17.

  Accordingly, such structural alteration should reduce the nucleophilicity of the serine residue, thus limiting the efficiency of the chemical process leading to the formation of the tetrahedral intermediate between enzyme and inhibitor. Ultimately, this drastic structural change reflects the steric clash that arises between the side chain of His468 and the restricted conformation of the alkyl chain imposed by the presence of the double bond between positions  and  for compound 15, which is nevertheless completely alleviated in 16.Even though compounds 14-16 can be viewed as very close analogues sharing similar intrinsic reactivity properties in the ketone moiety, the positional isomerism of the double bond has a dramatic impact on the inhibitory potency. According to present results, the apparently minor structural change associated with the presence of the double bond in position  or  is crucial, as it triggers a profound structural alteration in the catalytic triad upon formation of the tetrahedral intermediate with 15, while the adduct is properly accommodated upon positional isomerism of the double bond to position  in 16.

3 a

 3 IC 50 is the inhibitor concentration (M) that reduces the esterase activity by 50%. SE: standard error b This work; R = hexyl c R = butyl; taken from ref.[34]

Fig. 1 Fig. 2 Fig. 3 Fig. 4 13 Fig. 5 Fig. 6 Fig. 7 14 .Fig. 8 Fig. 9 Fig. 10

 123413567148910 Fig. 1Representation of selected ketone moieties in a series of TFMK inhibitors bearing a short alkyl chain (R is hydrogen in the model compounds used in quantum mechanical computations)

  Figure 1 Click here to download high resolution image

Figure 10

 10 Figure 10 Click here to download line figure: Fig10.doc

Table 1

 1 Free energy differences (kcal mol -1 ) for the formation of the hydrate species in the gas phase (G gas ) and in aqueous solution (G sol ) for model compounds 1-13

				Hydrate form
	Compound	G gas	G sol
				(%)
	1	-3.5	0.2	98
	2	-3.2	-0.2	99
	3	-3.6	-1.4	100
	4	-4.0	-0.5	99
	5	1.6	3.3	18
	6	-6.2	-2.8	100
	7	-3.1	-0.4	99
	8	4.5	4.7	2
	9	9.9	10.5	0
	10	2.0	2.5	44
	11	4.3	9.4	0
	12	2.9	6.3	0
	13	1.6	6.5	0

Table 2

 2 Experimental data for the inhibitory potency a of TFMK compounds 1-13

	Compound	-log IC 50 ( SE) b	-log IC 50 ( SE) c
	1	7.9 (0.04)	8.1
	2	7.4 (0.04)	6.5
	3	7.3 (0.04)	7.3
	4	7.7 (0.02)	7.3
	5	7.6 (0.06)	---
	6	5.9 (0.05)	---

Table 3

 3 Experimental data for the inhibitory potency (IC 50 ; nM) a of TFMK compounds 14-

	18
	Compound	IC 50 ( SE)
	14	14.0 ± 2.0
	15	> 1000
	16	11.8 ± 0.7
	17	9.9 ± 0.9
	18	> 1000
	a IC 50 is the inhibitor concentration (M)
	that reduces the esterase activity by
	50%. SE: standard error

Acknowledgments

We are indebted to Dr. H. Huang for technical assistance. We also gratefully acknowledge the Spanish Ministerio de Ciencia e Innovación (MICINN; projects AGL2006-13489-C02-01/AGR and SAF2008-05595) and the Generalitat de Catalunya (SCG2009 294 and 2009 SGR 871) for financial support. Partial support was provided by National Institute of Environmental Health Sciences R01 ES002710 and Superfund Research Program, P42 ES004699. JR acknowledges the Spanish Ministerio de Ciencia e Innovación for a FPI fellowship. BH is a George and Judy Marcus Senior Fellow of the American Asthma Association.

Supporting Information

Chemical characterization of compound 4 and its trifluoromethyl carbinol precursor. Plots of the time evolution of the potential energy and RMSD of the backbone atoms for the trajectories of the hCE1 complexes with inhibitors 14-16. * To whom correspondence should be sent: rpouplana@ub.edu, Phone: +34 934024557, Fax: +34 934035987 (RP), fjluque@ub.edu, Phone: +34 934024557, Fax: +34 934035987 (FJL)

Supporting Information

Reactivity versus

Supporting Information

Click here to download attachment to manuscript: SI_revised.doc

Chemical characterization of compound 4 and is trifluoromethyl carbinol precursor

Trifluoromethyl carbinol. 1 H NMR (300 MHz, CDCl 3 ):  4.11 (m), 3.70-3.42 (m, 4H), 3.17 (b, 1H), 1.58 (t, J = 6.8 Hz), 1.4-1.2 (m, 10H), 0.87 (t, J = 7.2 Hz) ppm. 13 C NMR (100 MHz, CDCl 3 ):  124.18 (q CF , J =280 Hz), 69.27 (q CCF , J = 30.6 Hz), 67.94, 31.77, 29.35, 29.32, 29.2, 25.91, 22.62, 14.04 ppm. 19