Jean-Baptiste Fiot

Laurent D Cohen

Parnesh Raniga

Jurgen Fripp

Efficicient brain lesion segmentation using multi-modality tissue-based feature selection and Support Vector Machines

Keywords: Image processing, Brain lesion, Segmentation, Classification, Support Vector Machines

Support Vector Machines (SVM) are a machine learning technique that have been used for segmentation and classification of medical images, including segmentation of white matter hyper-intensities (WMH). Current approaches using SVM for WMH segmentation extract features from the brain and classify these followed by complex post-processing steps to remove false positives. The method presented in this paper combines advanced pre-processing, tissue-based feature selection and SVM classification to obtain efficient and accurate WMH segmentation. Features from 125 patients, generated from up to four MR modalities (T1-w, T2-w, PD and FLAIR), differing neighbourhood sizes and the use of multi-scale features were compared. We found that although using all four modalities gave the best overall classification (average Dice scores of 0.54 ± 0.12, 0.72 ± 0.06 and 0.82 ± 0.06 respectively for small, moderate and severe lesion loads); this was not significantly different (p = 0.50) from using just T1-w and FLAIR sequences (Dice scores of 0.52 ± 0.13, 0.71 ± 0.08 and 0.81 ± 0.07). Furthermore, there was a negligible difference between using 5 × 5 × 5 and 3 × 3 × 3 features (p = 0.93). Finally, we show that careful consideration of features and pre-processing techniques not only saves storage space and computation time but also leads to more efficient classification which outperforms the one based on all features with post-processing.

INTRODUCTION

White matter hyper-intensities (WMH) are regions in the brain white matter (WM) that appear with bright signal on T2-weighted (T2-w) and fluid attenuated inversion recovery (FLAIR) MR modalities. They are a possible risk factor for Alzheimer's Disease (AD) and vascular dementia, with progression associated with vascular factors and cognitive decline [START_REF] Lao | Computer-assisted segmentation of white matter lesions in 3D MR images using support vector machine[END_REF]. To quantify these changes in large scale population studies, it is desirable to have fully automatic and accurate segmentation methods to avoid time-consuming, costly and non-reproducible manual segmentations. However, WMH segmentation using a single modality is challenging because their signal intensity range overlaps with that of normal tissue: in T1-weighted (T1-w) images, WMH have intensities similar to grey matter (GM), and in T2-w and proton-density (PD) images, WMH look similar to cerebrospinal fluid (CSF). FLAIR images have been shown to be most sensitive to WMH [START_REF] Anbeek | Probabilistic segmentation of white matter lesions in MR imaging[END_REF], but can also present hyper-intensity artifacts that can lead to false positives. To improve the WMH segmentation performance, additional discriminative information is extracted from multiple MR modalities.

The most successful lesion segmentation methods in the literature have been developed for the detection of multiple sclerosis lesions, with a recent grand challenge comparing the performance of various techniques [START_REF] Styner | 3D segmentation in the clinic: A grand challenge II: MS lesion segmentation[END_REF]. Lesion segmentation algorithms can be categorised into unsupervised clustering or (semi-)supervised voxel-wise classification [START_REF] Llad'o X | Segmentation of multiple sclerosis lesions in brain MRI: A review of automated approaches[END_REF]. Unsupervised methods suffer from the issue of model selection. Supervised methods such as neural networks [START_REF] Dyrby | Segmentation of age-related white matter changes in a clinical multi-center study[END_REF], k-NN [START_REF] Anbeek | Probabilistic segmentation of white matter lesions in MR imaging[END_REF], Naive Bayes classifier [START_REF] Scully | An automated method for segmenting white matter lesions through multi-level morphometric feature classification with application to lupus[END_REF] and Parzen windows [START_REF] Sajja | Unified approach for multiple sclerosis lesion segmentation on brain mri[END_REF][START_REF] Datta | Segmentation and quantification of black holes in multiple sclerosis[END_REF] have been proposed. Neural networks can be efficient but designing an appropriate network architecture and setting the parameters are difficult.

We present an SVM based segmentation scheme whose preliminary results were presented as a conference paper in [START_REF] Fiot | Efficient Lesion Segmentation using Support Vector Machines[END_REF], and inspired by the work in [START_REF] Lao | Computer-assisted segmentation of white matter lesions in 3D MR images using support vector machine[END_REF][START_REF] Zacharaki | Measuring brain lesion progression with a supervised tissue classification system[END_REF]. Lao et al. applied four steps: preprocessing (co-registration, skull-stripping, intensity normalisation and inhomogeneity correction), SVM training with Adaboost, segmentation and elimination of false positives. Our implementation utilises a similar but more advanced pre-processing pipeline and a simpler training procedure. As one of the primary causes of errors in other approaches is false positive cortical regions, we use information from multiple modalities to define a mask of potential WMH. This mask is built from patient specific tissue segmentation and atlas based population tissue priors. It leads to three main advantages compared to existing techniques. First, such careful feature selection enables to have a more accurate model without the use of boosting. Second it limits the areas where the classification is performed on the training set, which means a faster overall brain classification. Third, it reduces the false positive regions that are usually found with naive classifiers, so the advanced post processing required by other techniques [START_REF] Lao | Computer-assisted segmentation of white matter lesions in 3D MR images using support vector machine[END_REF] are not necessary. We also evaluated the relative value of each MRI acquisition protocol for WM lesion segmentation. This scheme is quantitatively validated on a significantly larger dataset including healthy aging, mild cognitive impairment and AD subjects. These results were compared with other supervised classification algorithms such as kNN, Naive Bayes, Parzen windows and decision tree.

CLASSIFICATION AND SUPPORT VECTOR MACHINE THEORY

Lesion segmentation can be formulated as a binary classification problem. The SVM technique [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond (Adaptive Computation and Machine Learning[END_REF] solves it in a supervised way: given l labelled features (x i , y i ) ∈ X × {-1, 1}, it builds a function f : X → R such that y(.) = sign(f (.)) is an optimal labeling function. The function f is a solution of the optimization problem:

f * = argmin f ∈H K 1 l l i=1 V (f (x i ), y i ) + γ f 2 K (1)
where K : X × X → R is a Mercer Kernel, H K its associated Reproducing Kernel Hilbert Space of functions X → R and its corresponding norm K , and V is the hinge loss defined as V (f (x), y) = max{0, 1 -y × f (x)}. The loss function V controls the labeling performance, and the second term controls the smoothness of the solution.

The optimization problem is convex because of the convexity of the hinge loss function. However as the objective function is not differentiable, the problem is reformulated with additional slack variables ξ 1 , . . . , ξ l ∈ R:

f * = argmin f ∈H K ξ1,...,ξn∈R 1 l l i=1 ξ i + γ f 2 K (2) subject to: ξ i ≥ V (f (x i ), y i ) ∀i ∈ {1, . . . , l}
The Riesz representation theorem states that the solution of (1) exists in H K , and can be written:

f * (.) = l i=1 α i K(., x i ) with α i ∈ R (3)
By plugging the expansion of f from (3) in (2), the optimisation problem becomes a finite dimension optimisation problem. Let the matrix K be defined as K i,j = K(x i , x j ). The optimisation problem is now: min α1,...,α l ∈R ξ1,...,ξ l ∈R 1 l l i=1 ξ i + γα T Kα subject to: (4)

ξ i -1 + y i l j=1 α j K(x i , x j ) ≥ 0 ∀i ∈ {1, . . . , l} ξ i ≥ 0 ∀i ∈ {1, . . . , l}
Let µ, ν ∈ R l be the Lagrangian multipliers. The Lagrangian of this problem is:

L(α, ξ, ν, µ) = 1 l l i=1 ξ i + γα T Kα - l i=1 µ i ξ i -1 + y i l j=1 α j K(x i , x j ) - l i=1 ν i ξ i (5) Solving ∇ α L = 0 leads to α * i (µ, ν) = yiµi 2γ ∀i ∈ {1, . . . , l}. Solving ∇ ξ L = 0 leads to µ i + ν i = 1 l .
The Lagrance dual function is:

q(µ, ν) = inf α,ξ∈R l L(α, ξ, ν, µ) = (6) l i=1 µ i -1 4γ l i,j=1 y i y j µ i µ j K(x i , x j ) if µ i + ν i = 1 l -∞ otherwise
The dual problem consists in maximising q(µ, ν) subject to µ ≥ 0, ν ≥ 0, and is equivalent to:

max 0≤µ≤ 1 l l i=1 µ i - 1 4γ l i,j=1 y i y j µ i µ j K(x i , x j ) (7)
Therefore the problem that α must solve is:

max α1,...,α l ∈R 2 l i=1 α i y i - l i,j=1 α i α j K(x i , x j ) = max α1,...,α l ∈R 2α T y -α T Kα
The training vectors with α i = 0 are called the support vectors. The optimization maximizes the margin, which is the distance between the decision boundary and the support vectors.

MATERIALS AND METHODS

Data

The dataset used in this paper comes from the AIBL study [START_REF] Ellis | The Australian imaging, biomarkers and lifestyle (AIBL) study of aging: methodology and baseline characteristics of 1112 individuals recruited for a longitudinal study of Alzheimer's disease[END_REF], where T1-w (160

× 240 × 256 image, spacing 1.2 × 1 × 1 mm in

Proposed algorithm

The proposed algorithm uses the standard supervised classification design for segmentation: given images and corresponding segmentations, the goal is to build a classifier to segment new images (Fig. 2). To obtain good performance, adequate preprocessing, mask and feature type have to be defined. This application-specific part is followed by a machine learning process.

As summarised in Fig. 3, the proposed algorithm consists of the following steps: Lesions can be seen in the FLAIR and T2-w as a bright signal.

Raw images

Manual segmentations

Input data

Pre-processed dataset

Application-specific

Mask definition

Feature definition

Training set

Machine learning

Classifier

Testing set Segmentation in the mask Preprocessing: Images were rigidly co-registered [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF], bias-field corrected [START_REF] Salvado | Method to correct intensity inhomogeneity in MR images for atherosclerosis characterization[END_REF], smoothed using anisotropic diffusion and histogram equalised to a reference subject. T1-w images were segmented into WM, GM, CSF using an Expectation-Maximisation approach with priors [START_REF] Acosta | Automated voxel-based 3D cortical thickness measurement in a combined Lagrangian-Eulerian PDE approach using partial volume maps[END_REF]. For each modality, features were extracted within the mask defined below, and scaled to [0 Mask creation: A global threshold on FLAIR images provides a high sensitivity, but poor specificity, which means it can only be used to define areas of interest. To further refine the areas of interest, we define the region W as the intersection of the dilated Colin WM mask [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF] (which was registered rigidly [START_REF] Ourselin | Reconstructing a 3D structure from serial histological sections[END_REF] then non-rigidly [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF] to the subject) and the WM mask (from the tissue segmentation in patient space). Using the mean µ W and standard deviation σ W of the FLAIR intensities on W , an intensity threshold of µ W + τ σ W on W (τ = 2 was the numerical value used) is used to define the mask M τ (Fig. ):

µ W = mean x∈W {F LAIR(x)} σ W = std x∈W {F LAIR(x)} I(x ∈ W ) = F LAIR(x)-µ W σ W Final Mask M τ Mτ = 1 if I(x) > τ 0 otherwise τ ∈ R
M τ (x) = 1 if F LAIR(x) > µ W + τ σ W and x ∈ W 0 otherwise
Machine learning: A subset of 10 000 features, with half belonging to the lesion class, the other half belonging to the non-lesion class, randomly selected and equally distributed among the training samples was used to generate the classifiers. A Matlab implementation solving SVM in its primal formulation was used [START_REF] Melacci | Laplacian Support Vector Machines Trained in the Primal[END_REF][START_REF] Melacci | Manifold regularization: Laplacian SVM[END_REF]. The chosen kernel was the (Gaussian) radial basis function. The width of the kernel and the regularisation weight were selected via a 10-fold cross validation. Then each image in the test set was segmented within the patient mask created. Pixels outside this region were set to the non-lesion class. As post-processing, all the connected components segmented as lesion with less than 10 voxels were removed. (with S the computed segmentation, GT the ground truth and λ counting the number of voxels in a volume), the number of true/false positive/negative (TP, FP, TN, FN) voxels, the specificity ( T N T N +F P ) and the sensitivity ( T P T P +F N ) computed on the full images. Higher is better for DSC, TP, TN, sensitivity and specificity. Lower is better for FP and FN. Statistical significance was analysed via the p-values of paired t-tests [START_REF] Ott | An Introduction to Statistical Methods and Data Analysis[END_REF]. We performed experiments to test the influence of the combination of modalities, the influence of the feature type and the influence of using the mask in pre-processing instead of in the post-processing. The performance of the SVM classifier was compared to the performance of other supervized classification algorithms. As the overall lesion load impacts the segmentation performance, as previously reported in [START_REF] Anbeek | Probabilistic segmentation of white matter lesions in MR imaging[END_REF], results are displayed for low (<3mL), moderate (3-10mL) and severe (>10mL) lesion loads. Finally, the impact of the parameter in the mask creation was evaluated by performance bounds of the segmentation performance.

RESULTS

Performance with regard to modality combinations

The segmentation performance was evaluated for various combinations of modalities (using 3 × 3 × 3 neighourhood features). Figure 5 shows DSC, FP and TP values for four single-modality and four multi-modality features. TN, FN, sensitivity and specificity values are similar for the various combinations, so corresponding graphs are not displayed. When using one modality, FLAIR gives the best performance. However, combining several modalities reduces FP and increases TP. Table I indicates that on low and moderate lesion loads the T1-w + FLAIR combination performs statistically better than FLAIR (T2-w + FLAIR does not). On the overall dataset, the T1-w + T2-w + FLAIR combination performs statistically better than FLAIR. The model with the four modalities performs the best (Fig. 5), but not significantly better than T1-w + FLAIR (p=0.50).

Performance with regard to features type

Using the four modalities, the segmentation performance was evaluated for various feature types (Fig. 6). With neighbourhood intensity features, a 5 × 5 × 5 size slightly increased the DSC compared to 3 × 3 × 3, but the difference was not statistically significant (p = 0.93). Pyramidal features with 4 dimensions did not perform as well as neighbourhood intensity features, but the DSC difference was not statistically significant (p=0.21 when compared with 3 × 3 × 3 features, p=0.18 with 5 × 5 × 5 features).

Performance with regard to the mask use

Using the FLAIR modality and 3 × 3 × 3 neighbourhood feature type, the impact of the mask on the segmentation performance has been evaluated. The use of the mask in the pre-processing instead of post-processing significantly decreased FP and led to a much better DSC (Fig. 7). The computation time in the prediction step being linear in the number of features to label, computing predictions for a significantly lower number of features (only within the mask) reduced the computation time (41 times computation speed-up on our dataset with τ = 2).

Performance comparison with other supervised classification algorithms

Using the FLAIR + T1 combination and 3 × 3 × 3 feature type, the performance of the SVM classifier was compared with several other supervised classifiers: kNN (with k=100 as in [START_REF] Anbeek | Probabilistic segmentation of white matter lesions in MR imaging[END_REF]), Naives Bayes, Parzen window and decision tree. The kNN classifier classifies a features x to the class with the highest cardinality among the k nearest neighbours of x in the training set.

ŷ(x) = argmax y Card{x i ∈ kN N (x) and y i = y} (8)
The Naive Bayes method computes the posterior probability of a feature x belonging to each class y, and classifies according the largest posterior probability (see Eq. ( 9)). To compute the parameters of the probability density of feature x given class y, the features are assumed conditionally independent given the class.

In the case of the Parzen window, the prior probability is estimated as

P (x|y) = 1 Card{i|y i = y} i,yi=y K(x, x i ) ( 10
)
where K is for example the Gaussian kernel. The Parzen window classifier is then found using the Bayes rule:

ŷ(x) = argmax y P (y|x) = argmax y i,yi=y K(x, x i ) i K(x, x i ) (11)
On this dataset, using combined 3 × 3 × 3 features from FLAIR and T1 modalities, SVM obtained the best DSC results, followed by Parzen windows, kNN, decision tree and Naive Bayes (Fig. 8). In terms of FP, kNN, Parzen window and SVM provide the best results. These tree classifiers also provide the best specifities. However, kNN and Parzen window classifier have a quite-low sensitivity, whereas SVM is most sensitive (followed by classification tree).

Performance bounds from mask parameter setting

The influence of the τ parameter in the mask creation on the segmentation performance was evaluated. Three performance bounds were computed, which are independent of the modality combination and feature type.

First, the upper bound of the final FP, which is equal to the number of FP for the "lesioneverywhere" classifier (worst-case scenario in terms of FP). Figure 9a shows the FP upper bound for different values of τ . The higher τ is, the smaller is the mask M τ (which is the area used for the machine learning), and therefore the lower is the FP upper bound. Second, the lower bound of the final FN, which is obtained when using the optimal classifier within the mask (best case scenario in terms of FN). Figure 9b shows the FN lower bound for different values of τ . The higher τ is, the smaller is the mask M τ , the more lesion features will be left out, and therefore the higher is the FN lower bound. Third, the upper bound of the final Dice score, which is similarly obtained using the optimal classifier within the mask. The higher τ is, the lower is the Dice upper bound (Fig. 9c). These graphs give insight on the impact of the τ parameter. If τ is too high, the final segmentation performance will be low, no matter how good is the classifier segmenting inside M τ . If τ is too low, the algorithm has a high risk of FP, which is a known drawback in WM lesion segmentation. Setting τ therefore involves a trade-off between the use of tissue-information to reduce the risk of FP and a high upper bound performance. The value τ = 2, which was selected for all experiments, is in the range of acceptable values decreasing the FP upper bound without decreasing too much the Dice upper bound.

CONCLUSION

We have presented a machine learning scheme applied to the WMH segmentation problem. Our approach is inspired by the previous work on SVM but has a number of differences. The value τ = 2 was the value selected for all experiments.

accurate segmentation results. Using our pipeline and our dataseet, SVM has a higher classification performance than other supervised algorithms such as kNN, Naive Bayes, Parzen window and decision tree.

In this work we also quantified the relative performance variations with regard to different modalities or feature types. Regarding the modalities, our results confirm that using all of the four modalities adds discriminative information and improves the segmentation results, as reported in [START_REF] Lao | Computer-assisted segmentation of white matter lesions in 3D MR images using support vector machine[END_REF]. However, our quantitative results show that using only FLAIR and T1-w can give similar performance at a lower cost. One reason could be the lower axial resolution of our T2-w and PD images. Regarding the features types, there is a trade off between the complexity, storage place and computation time versus the performance.

As other important contribution of this work, the mask we define and use in the pre-processing has several positive impacts. First, it improves the classifier performance as the training features are selected in regions of interest, which leads to better classifiers. We have given insight on the trade-off related to the threshold parameter selection in the mask creation. Increasing the threshold decreases the upper bound of FP (and therefore the potential risk of final FP). However, a low FP upper bound comes at a price as increasing the threshold increases the FN lower bound and decreases the Dice upper bound, which means a treshold too high would cause poor final performance no matter how good is the classifier. Second, computation time and storage space required are significantly lower (41 times lower on our dataset with the chosen threshold) as features and predictions are computed in a restricted area. Finally, using our mask in the pre-processing makes most of the complex postprocessing steps required in current state-of-art methods redundant. into the early detection and causation of AD. For a list of members of the AIBL Research Group, see www.aibl.csiro.au/partners.html.

  the sagittal, coronal and axial directions, TR = 2300 ms, TE = 2.98 ms, flip angle = 9 • ), FLAIR (176 × 240 × 256, 0.90 × 0.98 × 0.98 mm, TR = 6000 ms, TE = 421 ms, flip angle = 120 • , TI = 2100 ms), T2-w (228 × 256 × 48, 0.94 × 0.94 × 3, TR = 3000 ms, TE = 101 ms, flip angle = 150 • ) and PD (228 × 256 × 48, 0.94 × 0.94 × 3, TR = 3000 ms, TE = 11 ms, flip angle = 150 • ) were acquired for 125 subjects. WM lesions were manually segmented by one of the authors (PR), reviewed by a neuro-radiologist and used as ground truth in the classification (Fig. 1).

Figure 1 .

 1 Figure 1. Axial slices from one subject illustrating the different MR modalities and manual segmentation.Lesions can be seen in the FLAIR and T2-w as a bright signal.

Figure 2 .

 2 Figure2. Supervised classication algorithms for segmentation aim at building a classifier from images and corresponding segmentations. To obtain good performance, adequate pre-processing, mask and feature definitions have to be used. This application-specific part is followed by a machine learning process where a classifier is built from training examples, to segment new (testing) images.

Figure 3 .

 3 Figure 3. The proposed WMH segmentation pipeline is composed of 3 mains steps:pre-processing, mask creation and machine learning.

Figure 4 .

 4 Figure 4. The mask creation uses both FLAIR and T1 modalities to combine intensity-based and tissue-based properties. First, an Expectation-Maximisation (EM) technique on the T1 is used to generate WM/GM/CSF segmentation. On the intersection W of the patient WM and the registered Colin WM, a normalized scalar map is computed from the FLAIR intensities. A final threshold on this map provides the mask Mτ .

Copyright c 2012 Figure 5 .

 20125 Figure 5. Segmentation performance with different modality combinations (using the 3 × 3 × 3 neighbourhood intensity feature type).

Figure 6 .

 6 Figure 6. Segmentation performance with different feature types (using the 4 modalities).

Figure 7 .

 7 Figure 7. Using our mask Mτ in the pre-processing gives better results than using it only as a post-processing step.

Figure 8 .

 8 Figure 8. Segmentation performance with different algorithms (using the 3 × 3 × 3 neighbourhood intensity feature type).

Figure 9 .

 9 Figure 9. Performance bounds due to the threshold τ in the mask creation. Increasing τ decreases the upper bound of FP (a), increases the lower bound of FN (b), and decreases the lower bound of the Dice score (c).The value τ = 2 was the value selected for all experiments.

Table I .

 I p-values of paired t-tests using 3 × 3 × 3 features. Statistically significant differences (p < α = 0.05) in bold green.The dataset was randomly split equally into training and test sets. A classifier was built using the training set, and then used to segment the test set. Then training set and test set were swapped, another classifier was built, and the rest of the segmentations were computed. Results were then merged. Model performances were compared using the Dice score [20] DSC = 2 λ(S∩GT )

		Modalities	p-values of t-tests for lesion load in mL (number of subjects)
	Model 1	Model 2	< 3 (35)	3-10 (47)	>10 (43) Any (125)
	FLAIR	FLAIR, T2-w	0.47	0.19	0.62	0.38
	FLAIR	FLAIR, T1-w	0.047	0.032	0.57	0.070
	FLAIR	FLAIR, T1-w, T2-w	0.048	0.014	0.26	0.047
	FLAIR	FLAIR, T1-w, T2-w, PD	0.011	0.002	0.23	0.014
	FLAIR, T1-w FLAIR, T1-w, T2-w, PD	0.59	0.41	0.51	0.50
	3.3. Validation

λ(S)+λ(GT )

Copyright c 2012 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2012) Prepared using cnmauth.cls DOI: 10.1002/cnm

ACKNOWLEDGEMENT

We would like to thank Prof. Patricia M. Desmond from the University of Melbourne for her help with validation of the manual segmentations. Core funding for the AIBL study was provided by the CSIRO and was matched by contributions from the study partners (see www.aibl.csiro.au). The study also receives support from the National Health and Medical Research Council and the Dementia Collaborative Research Centres program (DCRC2). The Alzheimers Association (USA) has contributed support to allow AIBL neuroimaging stream data (and accompanying clinical and cognitive data) to be made available via the Alzheimers Disease Neuroimaging Initiative website. We thank Alzheimers Australia (Victoria and Western Australia) who have assisted with promotion of the study and screening of telephone calls from volunteers, and collaborated with AIBL volunteer functions. The AIBL team wishes to thank all those who took part as subjects in the study for their commitment and dedication to helping advance research Copyright c 2012 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2012) Prepared using cnmauth.cls DOI: 10.1002/cnm


Cristina Mazzoni
email: cristina.mazzoni@uniroma1.it

Eleonora Mangiapelo

Vanessa Palermo

Claudio Falcone

Manuela Côrte-Real

Michael Breitenbach

Paula Ludovico

Hypothesis: is yeast a clock model to study the onset of humans aging phenotypes?

Keywords: aging, apoptosis, cancer, colony, dedifferentiation, keratoses, yeast www.frontiersin.org

In this paper we report the growth and aging of yeast colonies derived from single cells isolated by micromanipulation and seeded one by one on separated plates to avoid growth interference by surrounding colonies. We named this procedure clonal life span, and it could represent a third way of studying aging together with the replicative life span and chronological life span. In this study we observed over time the formation of cell mass similar to the human "senile warts" (seborrheic keratoses), the skin lesions that often appear after 30 years of life and increase in number and size over the years. We observed that similar signs of aging appear in yeast colonies after about 27 days of growth and increase during aging. In this respect we hypothesize to use yeast as a clock to study the onset of human aging phenotypes.

Nowadays, Saccharomyces cerevisiae has been widely accepted as a model for the study of aging of multicellular eukaryotes [START_REF] Kaeberlein | Lessons on longevity from budding yeast[END_REF]. In fact, in this organism, we can measure the number of mitotic events an individual mother cell can undergo before senescence [START_REF] Mortimer | Life span of individual yeast cells[END_REF], referred as replicative life span (RLS). In addition, we can measure the time a non-dividing cell population can remain viable and this is called chronological life span (CLS; [START_REF] Fabrizio | The chronological life span of Saccharomyces cerevisiae[END_REF]. RLS has been suggested to be a model for the aging of mitotic tissues, whereas CLS has been likened to the aging of post-mitotic tissues [START_REF] Maclean | Chronological lifespan of stationary phase yeast cells; a model for investigating the factors that might influence the ageing of postmitotic tissues in higher organisms[END_REF] and both can induce apoptosis [START_REF] Buttner | Why yeast cells can undergo apoptosis: death in times of peace, love, and war[END_REF].

Different evolutionary studies are using yeast as a model to investigate on the initial emergence of multicellularity. The formation of multicellular aggregates in liquid cultures can be the result of incomplete cell separation or following the selection of clusters of cells whether by post-division adhesion or by aggregation [START_REF] Koschwanez | Sucrose utilization in budding yeast as a model for the origin of undifferentiated multicellularity[END_REF][START_REF] Ratcliff | Experimental evolution of multicellularity[END_REF]. On solid medium yeast forms colonies after repeated cell divisions.

Here we report the development and aging over 50 days of yeast colonies derived from individual cells of the wild type strain CML39-11A [START_REF] Mazzoni | Yeast caspase 1 links messenger RNA stability to apoptosis in yeast[END_REF] isolated by micromanipulation and placed each in one synthetic dextrose (SD) plate. We named this experimental approach clonal life span (ClLS) in that it allows to study the development of clonal cell avoiding the growth interference mediated by the acidic-base pulse generated by surrounding colonies [START_REF] Váchová | Aging and longevity of yeast colony populations: metabolic adaptation and differentiation[END_REF].

The organization of yeast colonies is ensured by signals transmitted and received by dividing cells within a colony and by chemical alkaline/acid pulse resulting from metabolic activity of colonies growing nearby.

Ammonia signaling is the first non-directed alkaline pulse produced by neighboring colonies and it is followed by a second step leading to acidification of the medium. A second ammonia pulse of higher intensity then occurs and is oriented toward the neighbor colonies. Ammonia signaling results in growth inhibition of the facing parts of near developing colonies [START_REF] Palková | Ammonia mediates communication between yeast colonies[END_REF].

The acid phase induces the production of reactive oxygen species (ROS) and other harmful products by all the cells forming the colony and induces apoptosis. The subsequent ammonia signal triggers metabolic changes that allow cells to lower their ROS production [START_REF] Palková | Life within a community: benefit to yeast long-term survival[END_REF].Thereishenceaselection of cells within the colony, determined also by the surrounding cells, which limit colony size and development. Cells located at the colony border grow slowly and are healthier compared to the ones located in the center of colonies, which predominantly undergo death.

In the absence of negative growth control exerted by neighboring colonies border cells of an individual colony can divide more and more times, probably without neutralizing ROS. Consequently, colonies might show larger size, peculiar phenotypes, and increased frequency of mutation during aging.

As shown in Figure 1, the size of the individually plated colonies increased within 50 days up to 12-13 mm. Interestingly, after 27 days we observed the appearance of some wart-shaped formations. The nature of these excrescences is still unknown, but one can hypothesize that during aging DNA damage and mutation frequency increase [START_REF] Fabrizio | Sir2 blocks extreme life-span extension[END_REF] leading to cells escaping growth control. Although not demonstrated for warts, a similar phenomenon occurs during yeast CLS where the regrowth of few cells, probably adapted mutants generated within aging populations, has been observed [START_REF] Fabrizio | Chronological aging-induced apoptosis in yeast[END_REF]. The latter authors suggested that such mutants are reminiscent of cancer cells, which become resistant to apoptosis and duplicate under conditions that are normally not permissive for growth. It has been demonstrated that these mutant cells arise from altruistic death program partially mediated by superoxide [START_REF] Fabrizio | Superoxide is a mediator of an altruistic aging program in Saccharomyces cerevisiae[END_REF]. The frequency of this cancer-like regrowth phenotype in yeast liquid cultures is greatly reduced under calorie restriction (CR) and in the presence of mutations in the Tor/Sch9 and Ras/AC/PKA (rat sarcoma/adenylate cyclase/protein kinase A) pathways. For these reasons yeast regrowth has been proposed as a useful phenomenon to study the age-dependent effect of mutations associated with cancer [START_REF] Madia | A simple model system for age-dependent DNA damage and cancer[END_REF].

Initially the yeast warts mainly appear in the middle of the colony, where older cells are located, and then they enlarge and propagate along the edge of the whole colony. Anyway, the number of such warts significantly differs from colony to colony suggesting the stochastic nature of these events.

To look closer at the nature of warts, we repeated the experiments extending the aging time up to around 100 days. We compared viability and mutation frequency of cells from warts (w) and from the smooth layer (L), picked from the same or different colonies. As shown in Figure 2A, w cells showed mainly higher viability compared to L cells. At the same time, we also determined the occurrence of mutations in both w and L cells by measuring the reversion frequency to the prototrophic phenotype of the auxotrophic mutation trp1-289 carried by the CML39-11A strain.

As shown in Figure 2C, the number of warts, compared to 50 days aged colonies, increased significantly, indicating that this phenomenon can occur for long time.

In Figure 2B is reported the number of revertants to Trp + phenotype (capability to grow in synthetic medium without tryptophan) normalized to 10 6 viable cells. As can be seen, the reversion frequency of w cells although heterogeneous, was higher than in L cells, reaching in some cases very high levels (i.e., w29, w31, and w32 values are out of scale).

In cancer cells, nuclear morphology is often altered and nuclei appear bigger and irregular in their contours [START_REF] Zink | Nuclear structure in cancer cells[END_REF]. 4 ′ ,6-Diamidino-2-phenylindole (DAPI) staining of DNA revealed the presence of round nuclei in cells coming from the smooth surface (Figure 3A) as well as in a fraction of w cells (not shown). In addition, part of w cells population showed abnormal nuclei morphology, including bigger dimension, fragmentation, and irregular contours (Figures 3B-D). Finally, the presence of many cells showing diffused DNA not organized in defined nuclei, together with multinucleated cells, indicates possible defects in cell division.

We like to speculate that these warts could represent the signs of age, just like senile warts [seborrheic keratoses (SKs)], the skin lesions that appear in humans around the age of 30 years and increase in number during aging [START_REF] Parish | Barnacles, old age marks, or just plain seborrheic keratoses[END_REF]. Seborrheic keratoses are the most common benign tumors in older individuals and they develop from the proliferation of epidermal cells. Although no specific etiologic factors have been identified, it is known that SK occurs more frequently in sunlight-exposed areas of the body such as the face and the neck [START_REF] Yeatman | The prevalence of seborrhoeic keratoses in an Australian population: does exposure to sunlight play a part in their frequency?[END_REF]. SKs are benign but secondary tumors, and Bowen disease (squamous cell carcinoma in situ) or malignant melanoma may occasionally arise within these lesions [START_REF] Terada | Pigmented Bowen disease arising in pigmented reticulated seborrheic keratosis[END_REF][START_REF] Böer-Auer | Cytokeratin 10negative nested pattern enables sure distinction of clonal seborrheic keratosis from pagetoid Bowen's disease[END_REF].

A similar situation seems to occur in yeast warts, as they can show both normal and cancer-like phenotypes.

Up to now, the patho-mechanisms of SK are not fully understood. Several studies showed that mutations in the fibroblast growth factor receptor 3 (FGFR3) and in the 110 kDa catalytic subunit of the phosphoinositide-3-kinase (PIK3CA) are present in human and mice benign skin tumors [START_REF] Logie | Activating mutations of the tyrosine kinase receptor FGFR3 are associated with benign skin tumors in mice and humans[END_REF]. More recently, some novel insight into the molecular basis of these benign skin lesions came from gene expression analysis study by DNA microarray that identified several upregulated genes, including the oncogenic form ΔNp63 of the transcriptional regulator p63 [START_REF] Seo | Microarray analysis reveals increased expression of DeltaNp63alpha in seborrhoeic keratosis[END_REF].

Recently, a yeast gene related to p63, NDT80, has been identified in controlling the aging process. In fact, this gene is involved in rejuvenilation of yeast cells during sporulation and yeast cells over-expressing this gene are able to double their lifespan [START_REF] Unal | Gametogenesis eliminates age-induced cellular damage and resets life span in yeast[END_REF].

We propose the ClLS approach as an additional method to study the onset of cells escaping growth control that, in turn, can be isolated from colonies for further studies, including genome wide analysis.

The treatment of SK human lesions varies from topical application of 5-fluorouracil (5-FU), cryotherapy, electrodessication, curettage to excisional surgery [START_REF] Park | Facial Plastic Surgery: The Essential Guide[END_REF][START_REF] Sand | Cutaneous lesions of the external ear[END_REF].

Yeast has already been successfully used to screen new compound and/or to assess their mechanisms of action [START_REF] Cassidy-Stone | Chemical inhibition of the mitochondrial division dynamin reveals its role in Bax/Bak-dependent mitochondrial outer membrane permeabilization[END_REF][START_REF] Regina | New arylthioindoles and related bioisosteres at the sulfur bridging group. 4. Synthesis, tubulin polymerization, cell growth inhibition, and molecular modeling studies[END_REF][START_REF] Palermo | Acetyl-Lcarnitine protects yeast cells from apoptosis and aging and inhibits mitochondrial fission[END_REF][START_REF] Palermo | Drug-induced inhibition of tubulin polymerization induces mitochondrion-mediated apoptosis in yeast[END_REF][START_REF] Palermo | Apple can act as anti-aging on yeast cells[END_REF]. In this respect, ClLS can provide a useful tool for the screening of molecules that are able to delay and/or reduce the onset of warts in yeast colonies.

In conclusions, we propose that yeast, by means of the regrowth phenotype during CLS and the appearance of warts during ClLS, can be a valuable model to study the formation of cell mass escaping the control of growth and it will be very interesting to study the nature of genes involved in this phenomenon. Moreover, yeast might represent a nice clock to measure the onset of aging phenotypes in humans.

ACKNOWLEDGMENT

This work was supported by PRIN 2009.

FIGURE 1 |

 1 FIGURE 1 | Growth and aging of colonies derived from of eight yeast cells of the wild type strain CML39-11A (Mat a, ade1-101, his3-Δ1, leu2, ura3, trp1-289) isolated by micromanipulation and placed each in one SD (0.67% yeast nitrogen base without amino acids), 2% glucose, and 20 µg/ml of the appropriate nutritional requirements according to the genotype of the strain) plate. Plates were incubated at 28 • C and recorded at the indicated days.

FIGURE 2 |

 2 FIGURE 2 | (A)Percentage viability of cells picked from single warts (w) or from the smooth layer (L) of colonies after 100 days of growth [shown in (C)]. About 1000 cells for each sample were analyzed for their capability to form microcolonies[START_REF] Palermo | Apoptosis and aging in mitochondrial morphology mutants of S. cerevisiae[END_REF]. The viability scale was fixed to 20% to better show the lower values. (B) Mutation frequency of the same cells measured as Trp + revertants normalized to 10 6 viable cells.

Frontiers in Oncology | Molecular and Cellular Oncology December 2012 | Volume 2 | Article 203 | 2

Frontiers in Oncology | Molecular and Cellular Oncology December 2012 | Volume 2 | Article 203 | 4

www.frontiersin.org December 2012 | Volume 2 | Article 203 | 5

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Received


Petr Popov

Sergei Grudinin
email: sergei.grudinin@inria.fr

Rapid determination of RMSDs corresponding to macromolecular rigid body motions

Keywords: RMSD calculation, rigid body motion, quaternion arithmetic, structure comparison, spatial transformations

Finding the root mean sum of squared deviations (RMSDs) between two coordinate vectors that correspond to the rigid body motion of a macromolecule is an important problem in structural bioinformatics, computational chemistry and molecular modeling. Standard algorithms compute the RMSD with time proportional to the number of atoms in the molecule. Here, we present RigidRMSD, a new algorithm that determines a set of RMSDs corresponding to a set of rigid body motions of a macromolecule in constant time with respect to the number of atoms in the molecule. Our algorithm is particularly useful for rigid body modeling applications such as rigid body docking, and also for high-throughput analysis of rigid body modeling and simulation results. We also introduce a constant-time rotation RMSD as a similarity measure for rigid molecules. A C++ implementation of our algorithm is available at http://nano-d.inrialpes.fr/software/RigidRMSD.

INTRODUCTION

The root mean square deviation (RMSD) is a widely used and powerful criterion to estimate the similarity between two ordered sets of points. In structural biology and bioinformatics, RMSD has been widely accepted as a measure of similarity between macromolecules. For rigid body modeling applications such as rigid body molecular docking 1,2 , rigid body molecular dynamics simulations 3,4 , and rigid body Monte Carlo simulations 5 , RMSD can be used as a measure of the rigid body motion of a molecule. However, determination of the RMSD can be a rate-limiting step for those applications where large number of rigid body motions should be compared. These applications range from conformation sampling in protein docking and structure-based drug design to high-throughput analysis of rigid body modeling and simulation results.

Much effort has been spent in developing algorithms for the optimal superposition of two molecules that minimizes the RMSD between the corresponding atoms [6][7][8][9][10][11][12][13][14][15][16][17][18] . In these methods, the squared RMSD is typically minimized with respect to the components of a rotation matrix or a rotation quaternion. However, in many applications of computational chemistry and structural bioinformatics a complementary problem emerges -given a set of rigid body motions of a reference molecule, compute the corresponding set of RMSDs. To the best of our knowledge, there exists no explicit description of an efficient algorithm for this problem in the literature. For the case of the RMSD between two positions of the same molecule after applying two spatial rigid body transformations, a formula can be found in the work of Rarey et al. 19 , however, it contains an error, which we correct below. Here, we present RigidRMSD, a new algorithm for constant-time RMSD computations. In particular, we provide a connection between the RMSD and the axis and the angle of the rotation. Also, we consider rotations represented by both matrices and quaternions, since the two representations are widely used in the description of spatial transformations. We demonstrate that the quaternion representation could be more efficient than the matrix representation. Our algorithm initializes in time linear in the number of atoms in the molecule and then computes the RMSD corresponding to a rigid body motion in constant time. Our algorithm can be very useful when computing multiple RMSDs corresponding to a sequence of rigid body motions, as e.g. in the DockTrina method 20 or clustering applications, as each new RMSD computation takes only constant time. To demonstrate the efficiency of the RigidRMSD library, we implemented an RMSD-based clustering algorithm and compared it with the standard clustering method. Finally, we provide several source-code examples that demonstrate the usage of our library.

THEORETICAL FOUNDATION

Weighted RMSD

Given a set of N points A = {a i } N and A ′ = {a ′ i } N with associated weights w = {w i } N , the weighted RMSD between them is given as

RMSD(A, A ′ ) 2 = 1 W i w i a i -a ′ i 2 , (1)
where W = i w i . Here, {w i } N are statistical weights that may emphasize the importance of a certain part of the structure, for example in case of a protein, the backbone or the side chains.

These weights can be also equal to atomic masses (in this case W equals to the total mass of the molecule) or may be set to 1 (in this case W = N ).

Quaternion arithmetic

A quaternion Q can be considered as a combination of a scalar s with a 3-component vector

q = {q x , q y , q z } T , Q = [s, q]. The product of two quaternions Q 1 = [s 1 , q 1 ] and Q 2 = [s 2 , q 2 ]
is a quaternion and can be expressed through a combination of scalar and vector products:

Q 1 • Q 2 ≡ [s 1 , q 1 ] • [s 2 , q 2 ] = [s 1 s 2 -(q 1 • q 2 ), s 1 q 2 + s 2 q 1 + (q 1 × q 2 )] . (2)
The squared norm of a quaternion Q is given as |Q| 2 = s 2 + q • q, and a unit quaternion is a quaternion with its norm equal to 1. An inverse quaternion Q -1 is given as

Q -1 = [s, -q]/ |Q| 2 .
A vector v can be treated as a quaternion with zero scalar component, v ≡ [0, v]. Then, a unit quaternion Q can be used to rotate vector v to a new position v ′ as follows

0, v ′ = Q [0, v] Q-1 = 0, (s 2 -q 2 )v + 2s(q × v) + 2(q • v)q = [0, v + 2q × (q × v + sv)] . (3)
Equivalently, the same rotation can be represented with a rotation matrix R, such that v ′ = Rv, where R can be expressed through the components of the quaternion Q as

R =    s 2 + q 2 x -q 2 y -q 2 z 2q x q y -2sq z 2q x q z + 2sq y 2q x q y + 2sq z s 2 -q 2 x + q 2 y -q 2 z 2q y q z -2sq x 2q x q z -2sq y 2q y q z + 2sq x s 2 -q 2 x -q 2 y + q 2 z    . (4)
A unit quaternion Q corresponding to a rotation by an angle α around a unit axis u is given as

Q = [cos α 2 , u sin α 2 ], and its inverse is Q-1 = [cos α 2 , -u sin α 2 ]
. Finally, N sequential rotations around different unit axes defined by unit quaternions

{ Qi } N result in a new vector v ′ according to 0, v ′ = QN QN-1 ... Q2 Q1 [0, v] Q-1 1 Q-1 2 ... Q-1 N -1 Q-1 N . (5)

Rigid body motion described with quaternions

Let R be a rotation matrix and T -a translation vector applied to a molecule with N atoms at positions A = {a i } N with a i = {x i , y i , z i } T , such that the new positions A ′ = {a ′ i } N are given as a ′ i = Ra i + T. Then, the weighted RMSD between A and A ′ is given as

RMSD 2 (A, A ′ ) = 1 W i w i |a i -Ra i -T| 2 . ( 6
)
We can rewrite the previous expression using quaternion representation of vectors a i and T as

RMSD 2 = 1 W i w i [0, a i ] -Q[0, a i ] Q-1 -[0, T] 2 . (7)
Here, the unit quaternion Q corresponds to the rotation matrix R. Since the norm of a quaternion does not change if we multiply it by a unit quaternion, we may right-multiply the kernel of the previous expression by Q to obtain

RMSD 2 = 1 W i w i [0, a i ] Q -Q[0, a i ] -[0, T] Q 2 . (8)
Using the scalar-vector representation of a quaternion, Q = [s, q], we rewrite the previous RMSD expression as

RMSD 2 = 1 W i w i [-q • T, -sT + (2a i -T) × q] 2 . ( 9
)
Performing scalar and vector products in Eq. ( 9), we obtain

RMSD 2 = 1 W i w i [q x T x + q y T y + q z T z ] 2 + [-sT x + q y (2z i -T z ) -q z (2y i -T y )] 2 (10)
+ [-sT y + q z (2x i -T x ) -q x (2z i -T z )] 2 + [-sT z + q x (2y i -T y ) -q y (2x i -T x )] 2 .
Grouping terms in Eq. ( 10) that depend on atomic positions together, we obtain

RMSD 2 = T 2 x + T 2 y + T 2 z + 4 W i w i q 2 x (y 2 i + z 2 i ) + q 2 y (x 2 i + z 2 i ) + q 2 z (x 2 i + y 2 i ) -2q x q y x i y i -2q x q z x i z i -2q y q z z i y i (11) + 4 W q x q z T z + q x q y T y -q 2 z T x -q 2 y T x + sq z T y -sq y T z i w i x i + 4 W q y q z T z + q x q y T x -q 2 x T y -q 2 z T y + sq x T z -sq z T x i w i y i + 4 W q y q z T y + q x q z T x -q 2 x T z -q 2 y T z + sq y T x -sq x T y i w i z i .
Introducing the inertia tensor I, the rotation matrix R, the center of mass vector C, and the 3 × 3

identity matrix E 3 , we may simplify the previous expression to

RMSD 2 = T 2 + 4 W q T Iq + 2T T (R -E 3 ) C, (12)
where C = 1 W { w i x i , w i y i , w i z i } T , rotation matrix R corresponds to the rotation with the unit quaternion Q according to Eq. ( 4), and the inertia tensor I is given as

I =    w i (y 2 i + z 2 i ) -w i x i y i -w i x i z i -w i x i y i w i (x 2 i + z 2 i ) -w i y i z i -w i x i z i -w i y i z i w i (x 2 i + y 2 i )    . (13)
Equation ( 12) is the principal result of this work. It consists of three parts, the pure translational contribution T 2 , the pure rotational contribution 4 W q T Iq, and the cross-term 2T T (R -E 3 ) C. In this equation, only two variables depend on the atomic positions {a i } N , the inertia tensor I, and the center of mass vector C. Below, we will use this fact when computing RMSDs for a set of rigid body motions.

RMSD corresponding to a pure rotation

An interesting consequence of Eq. ( 12) is the analytical expression of the RMSD for a pure rigid rotation. Recall that a unit quaternion in Eq. ( 12) can be represented as a rotation about a unit axis n by an angle α, Q = [cos α 2 , n sin α 2 ]. Then, if a rigid molecule is rotated about this axis passing through the origin, the RMSD for such a rotation is given as

RMSD 2 = 4 W sin 2 α 2 I(n), (14)
where I(n) is the reduction of the inertia tensor (13) to a scalar form for the unit axis n:

I(n) = n T In. (15)
The pure rotational contribution 4 W q T Iq in Eq. ( 12) can be rewritten in terms of a rotation matrix R as 4

W q T Iq = 4 W tr qq T I = 1 W tr (I) [1 -tr (R)] + 2 W tr (IR) . (16)
Here, rotation matrix R is connected with the vector part of the rotation quaternion q by Eq. ( 4).

Equivalently, Eq. ( 16) can be written as

4 W q T Iq = 2 W 3 i,j=1 (δ ij -R ij ) X ij , (17)
where matrix X is given as

X =    w i x 2 i w i x i y i w i x i z i w i x i y i w i y 2 i w i y i z i w i x i z i w i y i z i w i z 2 i    . (18)
Now, the weighted RMSD in Eq. ( 12) can be computed using the matrix description of the rotation:

RMSD 2 = T 2 + 2 W 3 i,j=1 (δ ij -R ij ) X ij + 2T T (R -E 3 ) C. (19)

RMSD corresponding to a relative rigid body motion

Let R 1 and R 2 be two rotation matrices and T 1 and T 2 -two translation vectors applied to a molecule with N atoms at positions A = {a i } N , such that new positions

A 1 = {a 1 i } N and A 2 = {a 2 i } N are given as a 1 i = R 1 a i + T 1 and a 2 i = R 2 a i + T 2 . Let a unit quaternion Q = [s, q] correspond to the relative rotation R T 2 R 1 .
Then, the weighted RMSD between positions A 1 and A 2 is given by a generalized version of Eq. ( 12) as

RMSD 2 (A 1 , A 2 ) = 4 W q T Iq + (T 1 -T 2 ) 2 + 2 (T 1 -T 2 ) T (R 1 -R 2 ) C. (20)
Using Eq. ( 17) we can rewrite the above equation using the matrix description of the rotation:

RMSD 2 (A 1 , A 2 ) = 2 W 3 i,j=1 δ ij - 3 k=1 R 1 ki R 2 kj X ij + (T 1 -T 2 ) 2 + 2 (T 1 -T 2 ) T (R 1 -R 2 ) C. (21)
The derived equation is equivalent to the formula obtained by

ALGORITHM IMPLEMENTATION Computational considerations

In the above equations (12-21), as we have mentioned above, only two variables depend on the atomic positions of the reference molecular structure -the inertia tensor I (or, its equivalent matrix X if the rotation is given by the matrix representation), and the center of mass vector Below we discuss computational strategies that allow to reduce the constant in O(N + M ).

In Eq. ( 12), the cross-term vanishes in the reference frame bound to the center of mass (COM)

of the molecule where C = 0. In this reference frame, the rotation is preserved, while the translation T COM is given as

T COM = RC + T -C. ( 22
)
We can equivalently obtain the translation in the COM reference frame using a rotation quaternion Q as

T COM = QC Q-1 + T -C. (23)
Therefore, in the COM reference frame, the RMSD can be computed with fewer arithmetic operations. More precisely, using the quaternion representation of the rotation, the RMSD is given as

RMSD 2 = T 2 COM + 4 W q T I COM q. ( 24
)
Similarly, if we use the matrix representation of the rotation, the RMSD is given as

RMSD 2 = T 2 COM + 2 W 3 i,j=1 (δ ij -R ij ) X COM ij . (25)
In the above equations, inertia tensor I COM and matrix X COM are computed in the COM coordinate system. A particularly interesting case is the computation of the RMSD in the principal axes of inertia (PAI) frame. The PAI frame is the coordinate system where the centre of mass vector C = 0 and the molecule is aligned along its principal axes, i.e. matrices I COM and X COM are diagonal.

In this frame, equations (24-25) are simpler. Also, in the PAI frame, RMSD corresponding to a relative rigid body motion defined by two rotation quaternions Q1 and Q2 and two translation vectors T 1 and T 2 will be

RMSD 2 (A 1 , A 2 ) = 4 W (s 1 q x 2 -q x 1 s 2 -q y 1 q z 2 + q z 1 q y 2 ) 2 I xx + (s 1 q y 2 -q y 1 s 2 -q z 1 q x 2 + q x 1 q z 2 ) 2 I yy + (s 1 q z 2 -q z 1 s 2 -q x 1 q y 2 + q y 1 q x 2 ) 2 I zz + (T 1 -T 2 ) 2 . ( 26
)
This equation uses three times fewer arithmetic operations compared to the previously published Eq. ( 21). More precisely, Eq. ( 26) requires only 38 arithmetic operations compared to 114 operations in Eq. ( 21). Generally, Eqs. (22-26) are more efficient in the number of arithmetic operations compared to Eqs. ( 12) and ( 21), as it is summarized in Table 1. This table lists the number of arithmetic operations needed to compute the squared RMSD using different representations of the rigid motion in three different coordinate systems, the world frame, the COM frame, and the PAI frame. As listed in Table 1, to compute the squared RMSD we need 54 arithmetic operations in the worst case, when the rigid rotation is given as a quaternion in the world frame. If we choose the coordinate system properly (the PAI frame), we can compute the squared RMSD in just 14 operations. Table 1 demonstrates that in the world frame one requires a fewer number of arithmetic operations to compute the RMSD if rotations are represented with rotation matrices, whereas in the COM and PAI frames the number of operations is equal between the two representations. However, when performing sequences of rotations, the quaternion representation is more numerically stable and computationally efficient compared to the matrix representation irrespective of the choice of coordinate system. Indeed, one requires 45 arithmetic operations to multiply two rotation matrices, whereas quaternion multiplication requires only 28 operations. Finally, Table 1 demonstrates that the squared RMSD for a relative rigid motion computed with the quaternion representation in the PAI frame requires three times fewer operations compared to the one computed with the matrix representation in the world frame.

Numerical tests

Throughout the article, we count the number of arithmetic operations in different equations according to the source code of the RigidRMSD library. We would like to mention that the cost of different arithmetic operations is not the same -division and square root are usually more expensive than multiplication, which is in turn more expensive than addition and subtraction [START_REF] Brent | Modern computer arithmetic[END_REF] . We should also mention that on modern computers minimizing the number of arithmetic operations is less important for the performance of a particular algorithm compared to increasing the amount of instruction level parallelism or improving memory access patterns and cache utilization, for example. Therefore, it is impossible to rigorously compare the performance of different algorithms solely based on their operation count. Thus, we only provide the total number of arithmetic operations as a rough estimation of the complexity of the equations and the corresponding algorithms. To get more practical numbers, in the following sections we run a series of tests with two different levels of compiler optimization.

We implemented the tests using the C++ programming language and compiled them using g++ compiler version 4.6 with optimization levels -O0 and -O3. For the gcc family of compilers, optimization option -O0 disables compiler optimization, whereas optimization option -O3 enables heavy optimization including interprocedural optimization and vectorization. We ran the tests on a 64-bit Linux Fedora operating system with Intel(R) Xeon(R) CPU X5650 @ 2.67GHz.

RESULTS AND DISCUSSION

This section presents numerical tests and practical applications of the equations derived in this article. First, we compare the quaternion representation with the matrix representation when computing sequential rotations (i.e. a composition of several rotations) and when computing a product of rotations with the subsequent RMSD computation. Second, we discuss the similarity measure between molecules and demonstrate that the rotation RMSD (Eq. 14) can be advantageous over a simpler angular distance measure. Finally, we present a rigid body clustering algorithm as an example of the application of the derived equations.

Rotation representation

Quaternions provide another way to represent rotations compared to conventional rotation matrices. In practice, the quaternion representation has several benefits over the matrix representation.

First, a quaternion compared to a matrix requires less storage, four values versus nine. Second, the orthonormalization of a quaternion costs much less than the orthogonalization of a matrix.

More precisely, orthonormalization of a quaternion can be accomplished by dividing the quaternion by its norm, which requires twelve arithmetic operations including one square root. However, there is no universal method for matrix orthonormalization. In this case, one may use the Gram-Schmidt orthonormalization method, QR decomposition, singular value decomposition or other methods, which are more computationally expensive compared to the quaternion orthonormalization [START_REF] Golub | Matrix Computations[END_REF] . Third, a product of two rotations using quaternions requires fewer arithmetic operations compared to the matrix representation (28 versus 45). Finally, the matrix multiplication is less numerically stable due to the accumulation of rounding errors. In summary, applications that require sequential rotations (e.g. some docking applications) will gain in speed, memory and numerical precision when using the quaternion representation.

To demonstrate the numerical efficiency of the quaternion representation, we ran a series of tests with two different levels of compiler optimization. In the first test, we performed 10 8 products of rotations using the two types of rotation representation and compared the timing for a single product of rotations with and without compiler optimization. The results of this test are presented in Table 2. We see that a rotation with quaternions is about 60% faster compared to a rotation with matrices regardless of the optimization level. In the second test, we computed a product of two rotations with the subsequent RMSD computation using Eqs. ( 22)-( 25) and repeated these operations 10 8 times. Then, we calculated the time required for a single product of rotations with the subsequent RMSD computation. The results of this test are also presented in Table 2.

Again, the quaternion representation is about 10% faster without optimization and 4% faster with optimization compared to the matrix representation. We should note that increasing the number of sequential rotations will provide a bigger speedup using the quaternion representation in this example. In the third test, we computed 10 8 RMSDs corresponding to a relative rigid body motion using the matrix representation of rotation (Eq. ( 21)) and the quaternion representation of rotation (Eq. ( 26)). We can see that our quaternion approach is 2.4-3.2 times faster compared to the matrix formula (Eq. ( 21)) depending on the level of compiler optimization.

To summarize, if a particular application operates with sequential rotations, as it happens in the DockTrina algorithm 20 or other docking applications, RMSD computations are more numerically efficient using the quaternion representation. Furthermore, the gain of using the quaternion representation is bigger up to 60 % when using a larger sequence of rotations.

Rotation RMSD as a similarity measure for molecular structures

It is still an open question how to measure the similarity between structures of a molecular complex [START_REF] Wallin | [END_REF] . For example, Rodrigues et al 24 developed a clustering method with the similarity measure based on the fraction of common contacts between two complexes. Another similarity measure was recently proposed by Vreven et al. 25 , where the angular distance computed in constanttime is used as the criterion for the similarity between the predictions from rigid body docking.

Nonetheless, the majority of the algorithms in the structural bioinformatics use the pair-wise RMSD as the similarity metric between the molecular structures.

Equation ( 14) is of particular interest when considered in relation to the aforementioned work of Vreven et al. 25 , where the authors demonstrated that the angular distance can serve as a similarity measure for rigid molecules as an alternative for the RMSD. More precisely, they defined the angular distance as the angle between the rotations corresponding to two docking predictions, ignoring the translational degrees of freedom. Vreven et al. claimed that the drawback of using the RMSD is that it is computationally expensive. However, we demonstrated that the RMSD can also be computed in constant time. Furthermore, in the context of Eq. ( 14), the angular distance is simply equal to the rotation angle α and does not take into account the geometry of a rigid molecule. In particular, for a fixed rotation angle, the angular distance for molecules of different size will be equal, while the RMSD can be very different. Another example that demonstrates the difference between the two measures is the rotation of a long linear molecule. The RMSD for such a rotation will dramatically depend on the axis of the rotation, while the angular distance will be the same regardless the rotation axis.

To conclude, we would like to emphasize that for comparison of rigid molecules of different size or molecules of non-spherical shape, it may be more rigorous to use the similarity measure defined by Eq. ( 14) instead of the angular distance. Particularly, our measure involves the scalar form of the inertia tensor (Eq. 15), thus taking into account the geometry and the rotation axis of the molecules.

Clustering algorithm

One of the possible applications of the RigidRMSD library can be the rigid body clustering. Molecular docking algorithms typically produce thousands of solutions, some of them having a very similar geometry. Therefore, it is practical to group these into clusters. As we have discussed above, there are multiple ways to measure the similarity between molecular structures [START_REF] Wallin | [END_REF] , however, most of the modern state-of-the-art clustering algorithms use the pair-wise RMSD as the similarity metric between the predictions, as it is implemented, e.g., in the Hex 1 and ZDOCK 2 docking algorithms. In the worst case, the complexity of such a clustering algorithm can be quadratic with respect to the number of docking predictions. Thus, an efficient pair-wise RMSD test can dramatically improve the clustering performance. The clustering algorithm used by the Hex and ZDOCK applications consists of the following steps. First, the docking prediction with the best score (yet unassigned to any cluster) is taken as the seed for the new cluster. Second, the pair-wise RMSDs between the seed and all other predictions (in case of ZDOCK) or some best predictions (in case of Hex) are measured and the predictions with the RMSD lower than a certain threshold are put into the cluster. Finally, these two steps are iterated until all docking predictions are assigned to corresponding clusters.

In order to demonstrate the efficiency of the RigidRMSD library, we compared the clustering algorithm implemented with our library to the one from the Hex software. We chose Hex for the comparison because it is a very fast rigid body docking tool and also because it explicitly provides the clustering time. It is worth to note that Hex's clustering algorithm has linear complexity with respect to the number of docking predictions, i.e. it is faster (though less accurate) than the standard RMSD-based clustering algorithms, as it is implemented in ZDOCK. Both Hex and ZDOCK clustering algorithms use the standard RMSD test linear in the number of atoms in the protein.

For the comparison, we collected a benchmark of 23 protein dimers of various size (see Table S1 in SI). After, we launched Hex version 6.3 on this benchmark and collected docking solutions before clustering, sizes of clusters, and clustering time. We then also clustered these solutions using the RigidRMSD library. Figure 1 shows the clustering time of the HEX clustering algorithm with respect to our clustering using Eqs. ( 21) and (26) as a function of the number of atoms in the smaller protein (left) and the number of docking solutions before the clustering (right). We can clearly see that our implementation of the clustering algorithm is more than an order of magnitude faster compared to the Hex implementation. Also, the quaternion representation of rotation, Eq.

(26), is on average three times more efficient compared to the matrix representation, Eq. ( 21). The efficiency of our clustering algorithm increases when using a larger RMSD threshold, as it is shown in Fig. S1 from SI. Also, mean cluster sizes obtained with our clustering algorithm are significantly larger compared to the Hex clustering (see Fig. S1 from SI), particularly at large RMSD thresholds.

This demonstrates that our implementation of the clustering algorithm is not only much faster, but also more accurate compared to the clustering in Hex, especially at large clustering thresholds.

CONCLUSIONS

We described a very fast and efficient way to compute the RMSD corresponding to the set of rigid body motions of a molecule. Our algorithm consists of an initialisation step followed by a series of constant-time RMSD computations. The initialization step has linear complexity with respect to the number of atoms in a molecule. However, each of the RMSD calculations requires only 14 to 54 arithmetic operations when using a single rigid motion (i.e. given with a single spatial rigid body transformation), or 38 to 114 arithmetic operations when using a relative rigid motion (i.e. given with a pair of spatial rigid body transformations), depending on the representation of the motion and the choice of the coordinate frame. This can be compared to 30 arithmetic operations needed to rotate a vector using a quaternion or 15 arithmetic operations needed to rotate a vector using a rotation matrix. We demonstrated that RMSD computations are more numerically efficient when using the quaternion representation of rotation. In particular, the gain of using the quaternion representation is bigger when using a larger sequence of rotations. We have also discussed two ways to measure the similarity between structures of a molecular complex.

In particular, we claim that it may be more rigorous to use the rotation RMSD similarity measure defined by Eq. ( 14) instead of the simpler measure based on the angular distance.

As an application of the RigidRMSD library, we implemented a clustering algorithm for solutions obtained with rigid body molecular docking tools. We showed that our implementation is more than one order of magnitude faster and also more accurate compared to the standard clustering algorithm used in the popular Hex docking software. A C++ implementation of our algorithm is available at http://nano-d.inrialpes.fr/software/RigidRMSD or by request from the authors. S1 in SI). For each protein complex, the number of considered docking solutions was fixed to 10,000. Right: Average time spent on clustering docking solutions by Hex and RigidRMSD with respect to the number of docking solutions. For this plot, we chose five structures with the number of atoms in the ligand protein of about 2,000 such that they result in a similar number of clusters and plotted the standard deviation of the clustering time for these structures. For both plots, time is plotted on a logarithmic scale and the clustering RMSD threshold is fixed to 10.0 Å. Table 2: Running time for three tests using two levels of compiler optimization. O0 optimization level disables optimization, whereas O3 optimization level enables heavy optimization including interprocedural optimization and vectorization. In the first test (columns 1 and 2), we performed 10 8 products of rotations using the two types of rotation representation and reported the timing for a single product of rotations. In the second test (columns 3 and 4), we computed a product of two rotations with the subsequent RMSD computation using Eqs. ( 22)-( 25) and repeated these operations 10 8 times for averaging. In the last test (columns 5 and 6), we computed 10 8 RMSDs corresponding to a relative rigid body motion, as in the clustering application, using the matrix representation of rotation (Eq. 21) and the quaternion representation of rotation (Eq. 26) and reported the timing for a single RMSD calculation.

Product

C.

  Therefore, given a set of M spatial transformations, we compute these two variables only once at the initialization step. The computational complexity of this step is linear with respect to the number of atoms N in the molecule. After, each RMSD computation for a single spatial transformation takes only constant time. The total cost to compute M RMSD values for a rigid molecule with N atoms thus will be O(N + M ), which is usually much smaller compared to the cost of standard algorithms, O(N M ), particularly at large values of M and N . More precisely, a standard algorithm computes the RMSD for each spatial transformation in O(N ) operations according to Eq. (1), thus resulting in O(N M ) overall complexity for M spatial transformations.

Figure 1 :

 1 Figure 1: Left: Time on clustering docking solutions by Hex and RigidRMSD with respect to the number of atoms in the ligand protein. Each point on the plot corresponds to a protein complex from the protein benchmark (see TableS1in SI). For each protein complex, the number of considered docking solutions was fixed to 10,000. Right: Average time spent on clustering docking solutions by Hex and RigidRMSD with respect to the number of docking solutions. For this plot, we chose five structures with the number of atoms in the ligand protein of about 2,000 such that they result in a similar number of clusters and plotted the standard deviation of the clustering time for these structures. For both plots, time is plotted on a logarithmic scale and the clustering RMSD threshold is fixed to 10.0 Å.

Table 1 :

 1 Number of arithmetic operations for the RMSD calculations with respect to different rotation representations and a different choice of the coordinate frame. These numbers were computed according to the source code of the RigidRMSD library. The references to the corresponding equations are given in the last column. These equations comprise only multiplication and addition / subtraction arithmetic operations.

		Multiplies Add/Subtract Total	Equation
	RMSD 2 (quaternion, world frame)	34	20	54	(24) and (23)
	RMSD 2 (matrix, world frame)	19	26	45	(25) and (22)
	RMSD 2 (quaternion, COM frame)	16	8	24	(24)
	RMSD 2 (matrix, COM frame)	10	14	24	(25)
	RMSD 2 (quaternion, PAI frame)	9	5	14	(24), I COM is diagonal
	RMSD 2 (matrix, PAI frame)	6	8	14	(25), X COM is diagonal
	RMSD 2 for clustering, (matrix, world frame)	55	59	114	(21)
	RMSD 2 for clustering, (quaternion, PAI frame)	21	17	38	(26)

ACKNOWLEDGMENTS

This work was supported by the Agence Nationale de la Recherche (ANR-2010-JCJC-0206-01 and ANR-11-MONU-006-01).


Fabien Todeschini
email: fabien.todeschini@st.com

Christophe Planat
email: christophe.planat@st.com

Philippe Bénabès
email: philippe.benabes@supelec.fr

EFFICIENCY THEORY OF SWITCHED CAPACITORS CONVERTER

Keywords: Wireless Sensor Network (WSN), energy harvesting, Maximum Power Point Tracking (MPPT), switched capacitor network

This paper reports on a theoretical study of the efficiency of a switched capacitor converter for autonomous sensors. The sensor is composed of an energy harvester (Solar cell), the electronic system (micro-controller, radio, temperature sensor, MEMS. . . ) and a battery.

Based on our theory, we calculate this efficiency with respect to: the switching frequency, the capacitance values and the battery voltage. We present this using 3-D graphs.

This theory is based on Ultra-Low-Power embedded systems such as wireless sensors with photovoltaic harvesters operating in low-light conditions. In this case, the harvested energy from PV do not exceed a few tens of micro-watts, which forbids the use of "power-hungry" coil-based power converters. On the contrary, using switched capacitors converter is a better alternative for this kind of application. We here present a theoretical framework deriving the maximum efficiency achieved by a switched capacitor converter with respect to: the switching frequency, the capacitance value and the battery voltage.

INTRODUCTION

Electronic autonomous devices are increasingly present in our everyday lives. These intelligent and smart systems can remotely collect data like temperature, humidity, sound, pressure, vibration, geomagnetic field. They are therefore convenient to acquire relevant physical knowledge about our surroundings. Because of the remote access constraint, wireless sensors need to harvest energy from their environment in order to extend their battery lifetime. The best way to extract the maximum amount of power from a harvester is to track the Maximum Power Point (MPP) [START_REF] Esram | Comparison of photovoltaic array maximum power point tracking techniques[END_REF]. Different techniques can be used depending on the required accuracy and/or the circuitry consumption and complexity constraints. The hill climbing and the Perturb and Observe methods are good examples to achieve a good accuracy but they need a microcontroller, which is not suitable for low power applications [START_REF] Esram | Comparison of photovoltaic array maximum power point tracking techniques[END_REF]. Faster techniques can also be used such as fractional open-circuit voltage [START_REF] Weddell | Ultra low-power photovoltaic mppt technique for indoor and outdoor wireless sensor nodes[END_REF]. The latter exploits the fact that the open circuit voltage (V OC ) has a near linear relationship with the voltage at maximum power point (V MP P ). These techniques require less processing but can result in inaccuracies as we do not know with certainty that we are working on the MPP.

Different kinds of power converters are used to extract energy depending on the harvested power and the harvester: the inductor based converters (buck, boost converters) and the switched capacitor converters (SC). Usually, inductor based converters are designed to convert high power (>1W) [START_REF] Ko Ko Win | Efficient solar energy harvester for wireless sensor nodes[END_REF] and require an external inductor. SC converters are more convenient as they can be entirely integrated on chip [START_REF] Seeman | An ultra-low-power power management ic for wireless sensor nodes[END_REF] and [START_REF] Shao | An inductor-less mppt design for light energy harvesting systems[END_REF].

But the efficiency of SC converters for Ultra-Low-Power remains low (below 50 % [START_REF] Maurath | A self-adaptive switchedcapacitor voltage converter with dynamic input load control for energy harvesting[END_REF]). In this paper, we propose how to calculate and evaluate accurately the achievable efficiency of SC converter charging a battery from a photovoltaic cell.

This paper is organized in the following manner: The charge pump principle is presented in section 2.1. Typical photovoltaic curves and models are discussed in section 2.2. The capacitor voltage is presented in section 2.3 and the efficiency calculation and the normalized efficiency are presented in section 2.4 and 2.5 respectively. The paper is ended by a conclusion in section 3.

SYSTEM DESCRIPTION

The charge pump principle

A charge pump consists of an n × m matrix of capacitors and some switches (Fig 1). The matrix topology changes at each clock cycle in order to adapt the harvester voltage level to the n . When this value is reached, the switches states change and the capacitors are discharged to the battery down to V batt m . The energy that is transferred from the PV to the battery is

E C = m × n × V batt m × ∆Q = V batt × C 0 V 1 - n m V batt (1)
The capacitor network acts as a voltage multiplier or divider, but unfortunately also acts as a resistor. The energy lost due to the capacitor discharge is proportional to

V 1 -nV batt m
. This loss essentially occurs when the PV voltage is much higher than the reported battery voltage n m V batt .

Figure 2 presents the capacitor voltage evolution. With a larger number of capacitors in the n × m matrix, the energy loss due to the capacitor discharge decreases. Indeed, because of the difference between the two voltages

VP V n -V batt m
if the value increases due to less capacitors (2 × 3 in this example), the loss will increase. On the contrary, with a lot of capacitors, the difference between V P V and V batt decreases, resulting in a smaller energy loss.

IV and PV modeling curves

In this example, the capacitor network is charged on the IV curve from the reported battery voltage n m × V batt that acts as a fixed voltage value, to the value V 1 . The voltage V 1 depends on the power (that is linked to the light condition), the capacitor value and the switching frequency. With a high switching frequency, V 1 will be close to n m × V batt . The same condition could be achieved with a smaller capacitor.

Thanks to these curves examples, we modeled them with a mathematical approach. This mathematical model is useful for the analysis and the matlab simulations especially because it is simplified. We modeled the IV curve in two parts. Indeed, the current does not change very much from the short circuit point to the optimum point. Therefore, we model this part by a small slope line starting from I SC . Then the curve decreases abruptly down to the Open Circuit Voltage. This is modeled by a negative parabola that decreases from the current at MPP to the V OC point. The Power vs Voltage curve is the IV curve multiplied by the voltage. The first part is the linear line starting from 0 up the Maximum Power Point. With the multiplication, the second part is a negative cubic curve up to the V OC . The Figure 4 shows the mathematical model of the IV and PV curves.

Capacitor voltage

The capacitors are charged while scanning the IV curve. The voltage law can be extracted from the differential equation

i = C du dt = IV (u × n) m ( 3
)
Where u is the voltage across a single capacitors, and IV is the current-voltage characteristic of the PV cell. Figure 5 represents the solution of this equation starting from 0V at t = 0. In reality, the capacitor's voltage evolves from V batt × n m to V 1 (being determined by the charge period).

First, as the current that flows through the capacitor is almost constant (i = I SC ) the voltage across the capacitor is quite linear:

time (t ISC C0 ) voltage (V) 0 5 10 15 20 1 2 3 n m V BAT T 5 V 1 V OC

Fig. 5:

Capacitor's voltage variation

u c = I SC + βt C t ≈ I SC C t (4)
Then the capacitors charge linearly. But as after the Maximum Power Point the current decrease abruptly down to 0, the slope of the voltage evolution across the capacitors decreases asymptotically toward V OC .

Efficiency calculation

The efficiency can be derived from the previous equations. The maximum energy that can be extracted from the P V during one period is at the M P P (we suppose that ∆T 2 in fig 2 is negligible compared to ∆T ):

E m = I mpp × U mpp × ∆T (5)
The effectively transferred energy depends on the value of the unitary capacitors, the sampling frequency (f s), and the battery voltage:

E t = C 0 × V batt × V 1 - n m V batt (6)
In this equation, V 1 is deduced from V batt and the switching period from the variation of the capacitor's voltage. (Fig 5)

The efficiency is then deduced as: η = Et Em .

We represent the evolution of the efficiency in a 3D plot (Fig 6). The X axis is the short circuit current divided by the product of the unitary capacitor value by the sampling frequency: Isc fs×C0 . The Y axis is the reported battery voltage normalized by U mpp .

The best efficiency (equal to 100 %) appears when the reported battery voltage is equal to the MPP and the frequency Fig. 6: Efficiency vs battery voltage and frequency x capacity is infinite. It means that this efficiency is theoretical and cannot be achieved in a circuit. This efficiency decreases rapidly if the frequency or the capacitor value decreases because we charge the capacitors with a value different from the Maximum Power Point. For the same reason, if we charge the capacitors from a value far from the optimum voltage, the efficiency will be low whatever the charging conditions (low/high frequency). In other words, if we start to charge the capacitors before the Maximum Power Point, reaching or exceeding this point is not useful and will not increase its efficiency.

Normalized efficiency

η ′ (V batt , f s) = η(V batt , f s) η(V batt , ∞) (7)
The next curve (Figure 7) shows the evolution of the normalized efficiency with respect to the frequency. We observe that when the battery voltage is low, the effi-ciency is not affected much by the frequency. It is then more efficient to have a low switching frequency if the starting voltage (V batt × n m ) is far from the Maximum Power Point. On the contrary, when the battery voltage is above the MPP, the efficiency decreases rapidly when the frequency decreases. Therefore, it is wise to avoid charging the capacitors after the Maximum Power Point.

CONCLUSION

These simulations have shown the theoretical efficiency that could be reached with a switched capacitor converter. The optimal switching frequency can easily be derived from these curves. In reality, the efficiency will be reduced due to the consumption of the switching logic that is directly proportional to the sampling frequency. Furthermore, designing a switched capacitor converter with a low switching frequency for low power applications is not interesting for high efficiency conversion.

Fig. 1 :

 1 Fig. 1: 2/3 switched capacitor principle

A

  solar cell is far to be a constant voltage source. To model the solar cell curves (Current vs Voltage and Power vs Volt-

Fig. 2 :

 2 Fig. 2: Capacitor voltage evolution of a 2x3 matrix

Figure 3

 3 shows the IV (Current-Voltage) and PV (Power-Voltage) characteristics measurements on a thin film silicon cell (amorphous solar cell). The IV curve starts from I SC and decrease slowly up to the voltage at Maximum Power Point (V MP P ). At this point, the curve drops abruptly up to 0A at the open circuit voltage (V OC ). The PV curve is the IV curve multiplied by the voltage. The Maximum Power Point (MPP) can be identified at the top of the PV curve. A thin film silicon solar cell has the advantage to present the same characteristics for indoor and outdoor light conditions (natural light like sun light or artificial light like incandescent bulb, fluorescent lights will not alter the shape of the curve). The I SC is proportional to the luminosity and the V OC level varies logarithmically with the luminance. With the thin film solar cell used for this work, the equation of the luminance is: luminance = 0.2819 ln(V OC ) + 2.7041

Fig. 3 :Fig. 4 :

 34 Fig. 3: IV and PV curve of an amorphous silicon solar cell

Fig. 7 :

 7 Fig. 7: Normalized efficiency vs battery voltage and frequency x capacity


F Todeschini
email: fabien.todeschini@st.com

C Planat

P Milazzo

S Tricomi

P Urard

P Benabes

A NANO QUIESCENT CURRENT POWER MANAGEMENT FOR AUTONOMOUS WIRELESS SENSOR NETWORK

Keywords: power management, wireless sensor network, energy harvesting, ultra-low power consumption

A Nano quiescent current Power Management for Autonomous Wireless Sensor Network dedicated to indoor and outdoor Photovoltaic (PV) Energy Harvesting is presented. The standby current consumption is 130nA @ 2.4V. The battery can be charged with less than 5% of the harvested energy from low light (70 lux) to high light (1 sun). The IP starts charging at light level as low as 10 lux providing ≃ 1µW . The voltage range is 2.4V -3.2V. The system load can be powered up to 400mA. This power management is manufactured using the CMOS 90nm technology in a 0.4mm 2 die area.

I. INTRODUCTION

Autonomous Wireless Sensor Network (AWSN) products are becoming possible thanks to the progress of ultralow power techniques in microelectronics leading to energy consumption reduction. One of the goal of AWSN is to avoid the use of disposable batteries, making it more eco-friendly. Another goal is to extend the system life time. In an AWSN node, the system load is composed of a micro-controller, one or several radios, and one or several sensors. Most of the time, this system consumes a greater current than the harvested one [START_REF] Zhi Ang Eu | Wireless sensor networks powered by ambient energy harvesting: An empirical characterization[END_REF]. So, the system has to be duty cycled and a rechargeable battery (or alternatively a supercap) has to be used as an energy tank. Thus, the battery is recharged by the harvester when the system is sleeping. The AWSN is viable when the energy levels are balanced between the quantity of harvested energy (including when the system is sleeping) and the energy consumed by the system (in active mode but also leakage and eventual counters inactive periods).

When the harvested energy is in the range of a few microwatts, the power management is not allowed to spend a large part of this energy. In this case, the use of hungry power converters is prohibited.

In [START_REF] Maurath | A self-adaptive switchedcapacitor voltage converter with dynamic input load control for energy harvesting[END_REF] a switched capacitor converter is used to convert the harvested energy to the load. This architecture can be fully integrated but no more than 50% of the input energy can be converted. In [START_REF] Yu | An ultralow input voltage power management circuit for indoor micro-light energy harvesting system[END_REF] a converter is needed to boost up the input voltage of the PV cell to charge the battery. These kinds of architectures are required when the harvester is working in a different voltage range than the battery. In [START_REF] Huang | A sub-100nA power management system for wireless structure health monitoring applications[END_REF] the circuit consumption is very low and the power management is used to power or standby the system from the battery when required. But this device is working only from battery to the system and doesn't manage or protect the battery. In [START_REF] Yifeng Qiu | 5µw-to-10mw input power range inductive boost converter for indoor photovoltaic energy harvesting with integrated maximum power point tracking algorithm[END_REF], an inductor based converter is used to convert photovoltaic energy harvested from indoor. The best end to end efficiency is 70% while consuming 0.8µA of static current. In [START_REF]Max17710 -energy harvesting charger and protector datasheet[END_REF], the charger uses a boost converter to charge a 4.2V battery but also requires a low-dropout regulator to power the system load if working at 3V. The use of 2 converters to power a load decreases the global efficiency. In [START_REF] Kadirvel | A 330nA energy-harvesting charger with battery management for solar and thermoelectric energy harvesting[END_REF] the power management is composed of an inductor based boost converter (DC/DC). With this converter, the input voltage cannot be equal or higher than the system voltage. This DC/DC powers directly the system and a battery or super-capacitor is used as an energy buffer. The buffer is protected against over and under voltage, but it's not protected against overcharging current. Some kinds of battery can be damaged with a high harvested current [START_REF] Maxell | Lithium manganese dioxide rechargeable battery[END_REF]. With this application, the voltage of the harvester should always be under the battery voltage. The reported power management quiescent current is at 330nA. A coil and other external components are required for system configuration.

In our proposal, the consumption of the power management is below 5% of the harvested current thanks to different working modes. These modes are autonomously managed. The three parts (the harvester, the battery and the system load) are wisely chosen to work in the same voltage range.

In this way, the efficiency can reach 87%. A Lithium-Ion battery charger, essential to protect the battery against high charging current, is also implemented.

This paper is presented in 5 sections: the circuit architecture will be presented in section II, the efficiency calculations are described in section III, the experimental results in section IV and the conclusion in section V.

II. POWER MANAGEMENT ARCHITECTURE

The power management (PM) has been designed to work with commercial batteries and PV. These parts are entirely part of the system. The photovoltaic cells are working from low light up to high light and are well suitable to work for indoor applications. We consider we will get a better efficiency without a hungry Maximum Power Point Tracking (MPPT) system.

Fig. 1 shows the block schematic of the power management. It can be supplied either by a 5V external power (i.e. USB) or by any source of power, in the range 0-6V. Depending on the light condition and the system current requirement, the PV can power the system and/or recharge the battery. This PM, supports batteries working in the range 2-3.2V. It embeds protection against over-voltage to avoid Li-Ion batteries over-charging. For this purpose, we use a Low Drop-Out (LDO) regulator that regulates the voltage at 3.2V. The under-voltage battery protection, designed to avoid the battery to be undercharged and be damaged, is designed with an hysteresis comparator that always senses the voltage of the battery. It will turns off the Power Switch (PS) when the battery is discharged to 2.4V. This will avoid the battery to be deeply discharge under 10% and then increase the battery life time (Lithium-Ion battery requirements). The comparator will turn PS when it is recharged to 2.8V.

We define different working modes to adapt to the light conditions and provide the best efficiency for the minimal power cost. The dark mode is the lowest power consumption mode. When no light is present, the system is powered only by the battery. The S2 switch, commanded by a resistor ladder that senses the PV voltage, is open to avoid leakage from battery through the PV LDO. The current flows from the battery through the Fast Switch (FS) that is switched on and through the switch PS. (S3 is consequently open). Only the battery protection and some very low power digital circuitry consume a part of the battery energy. The total consumption of the PM is 130nA at 2.4V.

The second mode is the low-light mode: from a few tens of lux to around 3klux. A voltage is present at the input of the PV LDO, the output of the ladder is higher than the threshold of the switch, S2 is then turned-on. The PV LDO is in follower mode, meaning that it is not regulating, the voltage drop is low (V OUT ≃ V IN ) and it provides all the energy harvested by the PV to the battery and/or the system. load current. The 2mA current limiter is off, saving the power and the battery is recharged by FS. The battery can also power the system through the switch PS. In this mode, the PV LDO consumes only 315nA at 150lux from the PV, which represent 1.75% of the harvested current. The total Power Management consumption is less than 580nA at 2.4V.

It also drives S3 (open) and FS (closed). This command occurs when the PV current capability is below the required

The third mode is the high-light mode. When the current from PV reaches a threshold fixed at 800µA the PV LDO switches in regulator mode, consuming 5µA, which represent 0.625% of the harvested current. The LDO regulates the internal common rail at 3.2V avoiding the battery to be overcharged. It switches on S3 and off FS. The 2mA limiter is active, consuming I charge 1000 . If the system requires a peak current higher than the PV can provide, FS switch is immediately closed in order to provide the required current. This is achieved by sensing V sys and the voltage of the battery. Once V sys < V batt (which means high current is required), the FS is active. It will be open when the system goes back to sleep. (V sys > V batt )

The last mode is the full-power mode. This mode occurs only when an external power is connected to the USB pad. As soon as a voltage is applied on it, S1 is switched on by the comparator and the USB LDO turns on and off the current limiter and the FS switch respectively and forces the PS switch whatever the battery voltage.

III. EFFICIENCY CALCULATION

The PV used for the design is working in the same voltage range than the battery. According to [START_REF] Belu | Design and analysis of a micro-solar power for wireless sensor networks[END_REF], the efficiency that we can have from a PV depending on the maximum power point and the effectively operating point is defined as:

η P V = P transf erred P MP P (1)
P MP P = V MP P × I MP P (2)
Depending on the current/voltage (IV) curves of the PV cell we are using, we can compute the power loss (or the efficiency) of the harvested power depending on the Maximum Power Point (MPP) and the Operating Point (OP), defined by V batt .

The figure 2 presents a typical measured current/voltage (IV) and power/voltage (PV) curve and the position of the MPP and OP on the Power vs Voltage (PV) curve. The table I presents the efficiency of the transferred energy for different light conditions at different battery voltages: The efficiency of this kind of power management is the ability to transfer the power. If the input power in entirely transferred to the output, the efficiency will be 100%. In other case, if there is a loss due to the joule effect, the efficiency will be:

voltage (V) current (µA) MPP OP power (µW) 0 1 2 V OP V MP P 4 V OC 0
η P M = P OUT P IN (3)
The efficiency of the power management has been computed for different values in the table II

IV. EXPERIMENTAL RESULTS

Simulations and chip measurements are matching and the power management has been tested in a global wireless sensor node. Table IV compares the power management consumptions between the simulations and measures for the tree harvesting modes: Table IV: Simulations vs Measurements.

❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵ ❵
The figure 3 shows the 3 modes of the power management: the dark mode, when the PV provides no current, the low-light mode where the current from PV goes directly through the battery and the high-light mode, the most power hungry mode, with regulation and current limiter for battery charging: The fig. 4 is the micrograph of the chip fabricated in 90nm with an effective area of 0.

I P V (

V. CONCLUSION

In this paper, a power management system from photovoltaic harvester to battery and to system load is presented.

Typical commercial power management doesn't include battery protection for charging current. Most of rechargeable lithium batteries are suitable for wireless sensors but require careful attention for charging. Our solution propose an areaefficient complete system power management for very lowenergy harvesting. With a suitable battery and PV, this PM is convenient for AWSN. The measurements show very low power consumption with less than 130nA standby current in a tiny 0.4mm 2 die area. The efficiency is higher than 85% from 30µA of the harvested current (corresponding to 200 lux: table II).

Fig. 1 :

 1 Fig. 1: Block schematic.

4 Fig. 2 :

 42 Fig. 2: IV and PV curve.

4mm 2 :Fig. 4 :

 24 Fig. 4: Micrograph of the Power Management in CMOS 90nm.

Table I :

 I Transferred energy on PV.

Table II :

 II PM efficiency.The global efficiency is calculated as the efficiency of the power management by the efficiency of the OP vs MPP:

	η = η P V × η P M		(4)
	η =	P transf erred × P OUT P MP P × P IN	(5)
	With this equation, we can compute the global efficiency
	from the previous tables. Different values are presented in
	the table III:
	P P P P P P lux V batt	2.4	2.6	2.8	3
	20	76.5%	72.5%	59.6%	33.7%
	200	77%	81%	83.5%	86.4%
	800	72.4%	77.7%	81%	85.3%
	2000	70%	74.6%	77.7%	82.8%

Table III :

 III Global efficiency.

Table V :

 V Comparison of this work with recently published work.