Face Recognition Using Invariance with a Single Training Sample

For the limits of memories and computing performance of current intelligent terminals, it is necessary to develop some strategies which can keep the balance of the accuracy and running time for face recognition. The purpose of the work in this paper is to find the invariant features of facial images and represent each subject with only one training sample for face recognition. We propose a two-layer hierarchical model, called invariance model, and its corresponding algorithms to keep the balance of accuracy, storage and running time. Especially, we take advantages of wavelet transformations and invariant moments to obtain the key features as well as reduce dimensions of feature data based on the cognitive rules of human brains. Furthermore, we improve usual pooling methods, e.g. max pooling and average pooling, and propose the weighted pooling method to reduce dimensions with no effect on accuracy, which let storage requirement and recognition time greatly decrease. The simulation results show that the proposed method does better than some typical and nearly-proposed algorithms in balancing the accuracy and running time.

Introduction
In recent years, face recognition technique has been applied on intelligent mobile terminals for varieties of applications. However, the terminals have some resource limits such as low power, small memories and low computing performance. In this case, the small training base and simple algorithms are preferred. So it is necessary to study strategies for face recognition with one training sample to keep a balance of the accuracy and running time. Actually, there are some fruits [1]- [5] in academics. The improved discriminative multi-manifold analysis (DMMA) method has been proposed in [1] by partitioning each face image into several non-overlapping patches to form an image set for each person. Although the recognition accuracy of DMMA [1] was higher than 80 percents for AR database, its running time on PC is as 100 times as typical methods such as PCA due to its complexity. Thus, one problem is addressed that whether an invariant factor we can acquire for face recognition using simple methods to perform less running time and an acceptable accuracy. This paper is trying to solve this problem.
The idea in this paper is inspired by the invariance in the visual recognition [6]- [8], which imitates the tuning properties of view-tuned cells in infero-temporal (IT) cortex [9], [10]. One of the corresponding famous models is Hierarchical Model and X (HMAX) model proposed by Poggio and his research group [11], [12]. HMAX was a four-layer hierarchical model, in which two types of computations: linear summation and non-linear max operation alternated between layers. HMAX model was first proposed for object recognition, not spatially for face recognition. The research group of Poggio later proposed the method for face recognition based on HMAX, which described cognitive characteristics in accuracy as a face descriptor [6]- [8]. This method is simulated on PC, not considering the operating platform, especially for the terminals without enough resources such as storages and computing abilities. The HMAX based methods mainly make best of Gabor filters for feature extraction, while Gabor filters could not get enough local details that represent the distinct features of different faces. Thus, we have to find a kind of wavelet filters instead of Gabor filters to get enough details on face features at different scales. To improve the accuracy using one training sample, the method of patch segments is also used for feature extraction. Furthermore, to reduce the quantities of data and still keep the invariance of the facial images, we take advantages of invariant moments [13], [14] and pooling techniques to keep the geometric invariance, such as rotation, shift and zooming, of the images as well as reducing the data storage.Based on these considerations, we propose an invariance model and its corresponding algorithms for face recognition using one training sample in conditions of low storage and low computing performance platform. This paper is organized as five sections. Section II introduces some related work including HMAX model and moment invariants. Section III illustrates the proposed invariance model and the corresponding algorithm in details. The simulation results based on public face databases are discussed in section IV. Section V makes a conclusion and discusses the further work.

Related Work 2.1. HMAX Model
There are two simple units and two complex units in HMAX model. The schematic of HMAX is shown in Figure 1. In Figure 1, S1 and S2 are simple units, and C1 and C2 are complex units. At S1, arrays of two-dimensional filters, generally Gabor filters, at four different orientations are used for the input raw images. Gabor filters are tightly tuned for both orientation and frequency but cover a wide range of spatial frequencies. Then, the groups of cells at same orientation but at a slightly different scales and positions are fed from S1 into C1 using MAX operation [11]. At S2,within each filter band, a square of four adjacent and non-overlapping C1 cells in 2×2 arrangement is grouped as cells, which are again fed from S2 into C2, finally achieving size invariance over all filter sizes in the four filter bands and position invariance over the whole input images.

Moment Invariants
Feature moment as a global invariant is often used for the feature selection to reduce the input for classifiers. The representations of seven moment invariants based on the secondorder and third-order normalized center moments are given by [7]. For a two-dimensional M× N image function f(x,y), the definitions of p+q order geometry moment mpq and the center moment are given by (1) and (2) .
(1) Here, is the gravity center on the direction of x in the image, is the gravity center on the direction of y in the image. Zero-order center moment is used to normalize the moments. So the normalized center moment is given by (3).
Then, the seven moment invariants are given by (4). The definition of the seven invariants are given in (4) for the second and third order moments. The seven invariants are useful for not only pattern identification independently of position, size and orientation but also independently of parallel projection. (4)

Invariance Model and Proposed Algorithms
As the illustration in section II, Gabor filters are generally used for feature extraction at the first layer in HMAX model. Gabor transformation could simulate the human visual system by decomposing the retinal image into a set of filter images. And the each filter image can represent changes in the intensities of the frequency and direction in the local scope of the raw image. So the texture features can be obtained by a group of the multi-channel Gabor filters. Equation (5) shows the Gabor filter, the key parameters are the frequency functions and the window size of Gauss function. Actually, Gabor transformation uses Gauss function as a window function to make local Fourier transformation by choosing the frequency and Gauss parameters. Although Gabor filters have tempo-spatial characteristics, they are not orthogonal so that different feature components have redundancies which lead to the low efficiency for texture feature extraction.
Here, sx and sy are the variances along x and y axes respectively, and f is the frequency of the sinusoidal function,  is the orientation of the Gabor filter.
We randomly select a facial image from ORL database shown in Figure 2.    Figure 3 shows that there are very slightly variances of the eight images and the more local details are not enough. On the other hand, for the same raw image, we make wavelet transformation at four scales using the wavelet 'db5'. Figure 4 shows the result, in which the upper left image, the upper right image, the lower left image and the lower right image are the detail images at scale four, scale three, scale two and scale one respectively. Comparing Figure 3 and Figure 4, wavelet transformation can obtain more details and local features than Gabor transformations. Actually, we let the most left image of the upper row in Figure 3 subtract the other seven images of Figure 3 respectively to get seven difference matrices, and then calculate the corresponding standard variance of the seven matrices. Similarly, we conduct the same operation on Figure 4. The computation results are shown in Figure 5. Obviously, the standard deviation values of the filtered images using wavelets are much bigger than that using Gabor. The standard deviation averages using wavelets and Gabor are 5.2420 and 0.7182 respectively. Therefore, we chose wavelet transformation instead of Gabor filters to get more local features. The typical orthogonal wavelet 'db5' is used as the kernel to make wavelet transformations. For each image, the wavelet transformation is conducted at four scales to get the corresponding approximation coefficients whose matrix is named A0 at the fourth scale and the detail coefficients whose matrices are named D1, D2, D3, and D4 at four different scales respectively. Then, A0, D1, D2, D3 and D4 are respectively singly used to make reconstruction so that to generate the reconstructed five images, e.g. Figure 4. After wavelet transformations and reconstructions, one image is decomposed to generate five images at different scales so that the quantities of data increased by five times. To solve this problem, the technique of dimension reduction has to be considered. Although PCA is a typical method of dimension reduction, it operates on the global features without considering the local details. While in case of only one training sample for each face, the local details play a very important role.
To keep local features as well as reduce data dimensions, we first divide the reconstructed images into patches. Combining the cognitive law and the facial image characteristics, we divide the facial image into 9 patches shown in Figure 6. Each patch just represents one physiological region of the face. Patch1 is a left forehead, patch2 is a middle forehead, patch3 is a right forehead, patch4 is a left cheek as well as a left eye, patch5 is a nose, patch6 is a right cheek as well as a right eye, patch7 is a left chin, patch8 is a lip, and patch9 is a right chin. Among so much feature data, the key features have to be preserved to represent the characteristics of each facial image. Because even if the images of the same face, they still have slightly variance in different conditions such as illumination, the invariance of features become more and more important. Thus, the technique of moment invariants is considered to preserve the invariance of feature data due to its invariance in scale, zoom and rotation. The details of moment invariants are introduced in section II.
For each reconstructed image, we first divide it into 9 patches, and then compute the seven moment invariants of each patch, and finally one column of 63 data is obtained. The corresponding algorithm is summarized in Table I. Step 2. Image reconstruction. We make wavelet reconstruction respectively using A0, D1, D2, D3 and D4 to generate five reconstructed images (e.g. Figure5).
Step 3. Patch segments For each reconstructed image, we divide it into 9 patches (e.g.Figure6).
Step 4. Moment invariants We calculate seven moments for each patch so that totally 9 groups of seven moments are obtained for each reconstructed image. Let the 9 groups of data be one column vector with the size of 63 data. Thus, one vector of 63 data represents the features of one reconstructed image (e.g.Figure7 and Figure8).
Step 5. Feature extraction According to the results of step 4, there are 63 feature data for one reconstructed image. Because one facial image is corresponding to five reconstructed images, there are totally 315 feature data for one facial image. We randomly chose some images shown in Figure 7 and Figure 8 to illustrate feature extraction. Figure 7 shows that the moment invariants of the two faces of the same person have only slight differences, especially, among the seven moment invariants referring to equation (4), 1  and 2  are almost the same. While Figure 8 shows the moment invariants of the two different facial images. Obviously, there are differences between the circle lines and the cross lines, especially in the reconstructed images at scale 1, 2 and 3(see Figure 8 To reduce the dimensions of feature data furtherly, we also proposed a weighted average pooling method. Observing the patches in Figure 6 and combing the cognitive rules of human brains, we suppose patch 4, 5, 6 and 8 have more contributions than other patches for representing one face, so we distribute comparatively bigger weights named 'w4','w5','w6','w8' to patch 4, 5, 6 and 8 than those named 'w1','w2','w3','w7','w9' to patch 1,2,3,7,9. All of the nine weights are experiencing values from many experiments, and the sum of all the weights is equal to 1. The details are shown in Table II. Using the weighted average pooling method, each reconstructed image has seven feature data so that each train sample has 35 feature data. The dimensions of training data are greatly reduced. That is, the storage requirement decreases which is a very important advantage for embedded device, e.g. sensor nodes of small size and low cost. Step 1: Weights distribution Distributing 'w 1 ','w 2 ''w 3 ','w 4 ','w 5 ','w 6 ','w 7 ','w 8 ','w 9 ' to patch 1,2,3,4,5,6,7,8,9 respectiv sum of the nine weights is 1.
Step2 :Weighted average pooling Supposing the moment invariants vectors of the nine patches are respectively v 1 ,v 2 ,v 3 ,v 4 ,v 5 ,v 6 ,v 7 ,v 8 , and v 9 , in sequence, then the feature vector of each reconstr is In a summary, an invariance model of two layers is constructed. Figure9 shows the invariance model structure composed of invariance attributes (IA) layer and invariance cluster (IC) layer. In IA, wavelet transformations are conducted at four scales resulting into four-scale detail coefficients and approximation coefficients, which are respectively used to generate new five reconstructed images. To keep the invariance and reduce quantities of data, we take advantages of the technique of patch segments and invariant moments to keep the global invariance of each reconstructed images. The corresponding algorithm of AI is summarized in A.1. On the base of AI, feature cluster has to be conducted to reduce dimensions as well as keep feature invariance. So in IC, the improved weighted pooling technique is used to select the key features. The corresponding algorithm of IC is summarized in A.2. Finally, a vector of 35 feature data represents the raw image.

Simulations
Three public face databases: ORL database, AR database and FERET database are used for simulations [15]. The details of selected facial images are shown in Table II. The simulation is running on PC with 3.3GHz CPU and 4G RAM. We first chose the proposed method, called IM short for invariance model, including A.1 and A.2, to make feature extraction and dimension reduction, and then chose KNN and LDA, the simple classifiers to make classification. To evaluate the proposed method, after step 2 in A.1, we use PCA, the typical dimension reduction technique to make dimension reduction instead of moment invariants, and then chose the same classifiers to make classification.      Figure 10 shows the recognition rate of the proposed method using classifiers of KNN and LDA respectively. Obviously, KNN and LDA have almost the same recognition rate for the three databases. Especially, the proposed method achieves the best accuracy of about 85% for AR database. This is because there are all frontal faces in AR database while in ORL and FERET databases there are some slight tilt faces. Figure11 and Figure12 show that using the same classifier, the recognition rate of the proposed method is much better than that of PCA in ORL and FERET databases. While in AR databases, the accuracies of both methods are almost the same. This also estimates that the proposed method is more robust than PCA in case of non-frontal faces. Similarly in Figure 13, KNN and LDA have slight effect on the accuracy for the both methods.
To evaluate the recognition time of each facial image, we also calculate the recognition time for each test image supposing the training database is ready. Table IV shows the average recognition time of each test image in the three databases. Obviously, the proposed method running time is generally less than PCA.

Conclusion
In this paper, we discuss the proposed invariance model. In this model, wavelet transformations are used instead of Gabor filters by analyzing the reconstructed images from decomposition coefficients to get more accurate local details. To avoid the effects caused by shift, rotation and other inference, the moments and weighted pooling technique are used to keep the invariance and reduce dimensions greatly. The simulation results show the proposed method keeps the balance of the accuracy, storage and running time.


extfieldsofstudy:
Biology
Computer Science

text:
Numerical modeling of inhomogeneous DNA replication kinetics

We present a calculation technique for modeling inhomogeneous DNA replication kinetics, where replication factors such as initiation rates or fork speeds can change with both position and time. We can use our model to simulate data sets obtained by molecular combing, a widely used experimental technique for probing replication. We can also infer information about the replication program by fitting our model to experimental data sets and also test the efficacy of planned experiments by fitting our model to simulated data sets. We consider asynchronous data sets and illustrate how a lack of synchrony affects replication profiles. In addition to combing data, our technique is well-adapted to microarray-based studies of replication.

INTRODUCTION
The development of epifluorescence microscopy techniques in combination with single-molecule methods such as molecular combing has had an enormous impact on the study of DNA replication [see Refs. [1][2][3] for reviews]. Earlier approaches such as autoradiography [4] and 2D gel electrophoresis [5] were not only slow but also provided only limited and localized information about the replication process. The advent of reliable, highthroughput methods such as molecular combing led to genome-wide studies of DNA replication kinetics. Moreover, the possibility of genome-level observations made possible complex experiments such as a study of the checkpoint regulation of DNA replication [6] and one of the relation between replication fork movement and origin spacing [7].
Molecular-combing measurements provide information about the DNA replication kinetics at a given time during the S phase [8]. These snapshots of the replication process can be used to improve our knowledge about the location of replication origins, along with their activation time [6,7,[9][10][11][12][13][14][15][16][17][18][19][20]. Such an analysis of the origin-firing program can, at least qualitatively, be based on a simple interpretation of combing measurements. For example, even when studying cell cycles that are completely asynchronous, we can infer from a series of combing measurements the relative replication times of different regions of the molecule. This is of biological importance when considering both local and global changes in replication patterns that occur during normal cell differentiation, oncogenesis and oncogene-induced senescence [21,22], Pluripotent embryonic stem cells, for example, rely on smaller replicon sizes that allow for diversification into distinct somatic cell lineages [23]. Knowledge of how replication patterns change during differentiation can give insight about underlying regulatory mechanisms and have contributed to our understanding of the impact of chromatin structure on technologies such as cloning [24]. Experiments based on synchronized cell populations can shed even more light on the origin positions and firing times, providing a detailed description of the space-time initiation program of DNA replication.
In order to fully exploit the large data sets obtained from molecular-combing experiments, we must rely on computer-modeling techniques. An approach commonly used to interpret combing results is to perform Monte Carlo simulations of hypothetical replication scenarios [25][26][27][28][29]. In such cases, the models are not designed to reproduce all the biological details of DNA replication but rather are focused on the kinetic aspects of the process (i.e., the rate at which origins are activated or the speed at which DNA is duplicated by replication forks). Nevertheless, fundamental biological features relating to the effect of chromatin modulations on replication patterns can be readily investigated. Therefore, these simulations allow one to test hypotheses about the overall replication scenario and to quantify the kinetic parameters of a given scenario. However, such simulations can be computationally intensive and therefore too slow to explore a series of scenarios (e.g., when performing a fit). The numerical modeling approach we present here addresses this issue by allowing us to study the average replication kinetics without having to do simulations. More precisely, our calculation technique gives the average replication rates and forks densities one would measure if an infinite number of replication cycles were observed or simulated. The speed gain from using our method compared to simulations makes it possible to explore a variety of kinetic parameters and to fit experimental data. This will facilitate automation of sample analysis that is of potential clinical and diagnostic interest because many genetic diseases such as cancer involve altered cell cycle and DNA replication kinetics [30].
In the following sections, we will present our calculation technique and show how it can be used to study a given replication scenario. Our modeling approach is flexible enough to accommodate inhomogeneous DNA replication kinetics, which are scenarios in which the initiation rate and replication fork speed can change with position (and time) along the molecule. Indeed, the methods presented here were motivated by single-molecule data that show clearly that the initiation of DNA replication in mammalian cells can occur in zones and that these zones are under developmental control [11]. In addition, large-scale spatial variation of initiation rates and developmental control of such rates has been seen in many microarray-based experiments that give information about local, population-averaged replication rates [for example, [31]]. Different chromatin environments, for example, correspond to regions with significantly different replication kinetics. The type of analysis discussed here can facilitate experiments addressing questions as to which chromatin factors (acetylation, methylation, phosphorylation) are involved in determining, on a concept specific basis, fork rates and origin efficiencies over extensive genomic regions.
In the Materials Section, we will describe the type of combing data considered in the examples presented in this article. The Methods Section will present the details of our calculation, as well as a comparison with simulated data. This section will also show how our model can be adapted to study asynchronously replicating cells and how we can infer information about the replication program by fitting simulated combing data. Finally, in the Notes Section, we will briefly discuss how a global analysis of all the replication kinetic parameters could be performed. We will also discuss how the effects of fork blocks (or pause sites) can be included into our model. This is also of potential clinical interest when examining the effects of a drug on DNA replication kinetics. Recently, for example, it has been shown that sister replisomes decouple when they encounter damaged DNA [32]. Therefore, the number of unidirectional forks might be an indicator of the number of blocked replication forks in a normal cell cycle. Such forks can be used as a marker to quantify and monitor transient DNA damage and DNA repair efficiencies on a genome-wide basis [33]. This is just one example among many potential applications of this type of modeling.

Kinetic description of DNA replication
The experimental situation we consider in this article is sketched in Fig. 1. Part (a) is a space-time illustration of the DNA replication kinetics during one cell cycle. In this representation, replication can be initiated at various positions x along the genome and at different times t since the beginning of S phase (see the circles in Fig. 1(a)). A pair of oppositely moving replication forks emerges from each initiation site, and each fork starts to duplicate DNA at a speed v(x, t). The dependence on x and t means that v can change with both the location along the genome and the time during S phase. For simplicity, the fork velocity used in this article is assumed constant throughout both space and time, as illustrated in Fig. 1. The forks progress until they reach one end of the molecule or until they coalesce with an oppositely moving fork (diamond symbol in Fig. 1(a)).

Molecular labeling
The labeling procedure consists in supplementing the regular nucleotides used for the DNA synthesis with nucleotide analogs that can be detected by immunofluorescence microscopy. If the cell cycles are synchronized, the nucleotide analogs can be introduced at a precise time t i after the beginning of S phase. Figure 1(b) shows examples of fluorescently labeled molecules that might be collected given the replication cycle presented in Fig. 1(a). Each example, from bottom to top, illustrates the labeled molecule obtained if the analogs were added at times t 0 , t 1 , t 2 , t 3 , t 4 , t 5 or t 6 . Basically, this procedure labels the DNA that is replicated after the specific time t i . In an actual experiment, a second type of nucleotide analog would be used to map the rest of the molecule and hence obtain a fully labeled molecule. In the simulations presented in this article, a second label is unnecessary, since there is no ambiguity about the molecule alignment or stretching. The results we will present here focus on the replication fraction of a molecule, f (x, t), which is the probability that the section of the molecule located at the position x along the genome is replicated by a given time t during S phase.

Data collection
In order to calculate the replication fraction f (x, t), we must ensure that the experimentally collected labeled molecules meet the following conditions: • the precise time when nucleotides were added must be known; • the molecules must be combed and stretched to identify the location of the labeled regions; • the molecules must be oriented (which can be achieved using hybridization probes).
Once these requirements are fulfilled, as in Fig. 1(b), we can compute the replication fraction profile f (x) of each labeled molecule by assigning the value f = 1 (replicated at time t) to the non-labeled regions of the molecule and f = 0 (unreplicated at t) to the labeled ones. Figure 1(c) presents the complete space-time replication profile resulting from our labeling procedure for the replication cycle example shown in Fig. 1(a). The labeling sequence can also be used to locate the positions and directions of replication forks as a function of time. Although it is possible to analyze fork data, too, we focus here on the replication fraction [see Ref. [34] for more on analyzing fork-density data].
In the following sections, we will compare calculated replication fractions to simulated ones. Our simulated data sets mimic the data one could experimentally collect, in the sense that we simulate each step of a similar combing experiment. On the other hand, our simulated data sets are "better" in some ways than their experimental counterparts: • the simulated cell cycles are perfectly synchronized; • the time when nucleotide analogs are introduced is known exactly (and once added, labels are instantly used for sequencing); • the positions along the simulated molecules are clearly defined (simulated results represent molecules that are perfectly stretched, aligned, and oriented).
Such ideal simulated data allow us to obtain precise statistics about the kinetics of replication itself, independent of experimental noise. (We could also simulate the noise processes, but that would complicate the presentation here.) As we will see below, these simulations can be used to test our numerical calculations.

Monte Carlo simulations of combing experiments
As mentioned above, we need to generate simulated data in order to test our numerical model of DNA replication. The simulations we carried out follow the requirements listed in Section . For each simulated labeled molecule used in this article, we performed a Monte Carlo simulation of the replication dynamics as presented in Fig. 1(a). The numerical-simulation procedure we followed is presented in detail in Ref. [35] and can be summarized as follows: 1. Define a simulation time step δt.
2. Begin a cell-cycle simulation with an unreplicated molecule of size L at time t = 0, which marks the beginning of S phase.
3. Simulate the initiation of replication origins (circles in Fig. 1(a)) between t and t + δt using a predefined initiation rate I(x, t). This rate is defined as the number of initiations per time per length of unreplicated DNA.
4. Allocate to each origin that fires a pair of oppositely moving replication forks (triangles in Fig. 1(a)).

5.
Calculate the displacement of all existing replication forks during the time step δt using the fork velocity profile v(x, t). Fork movement corresponds to DNA replication. 6. Delete forks that coalesce with oppositely moving forks (diamond in Fig. 1(a)) or that reach the end of the molecule.

Repeat from
Step 3 until the replication of the molecule is completed.
8. Define the time when the fluorescent nucleotides were introduced and deduce the corresponding labeled molecule pattern (as illustrated in Fig. 1(b) for various label addition times). As discussed below, this step is the one when we distinguish between synchronously and asynchronously replicating cells.
9. Repeat from Step 2 to collect more molecules.
The simulation steps described above are valid if all cell cycles are perfectly synchronous. For asynchronously replicating cells, we need to estimate the distribution of S phase start times (for example, using a Gaussian probability function). Then, cell cycles themselves can be simulated following Steps 1-7 above. However, the time when fluorescent nucleotides are added must be corrected relative to the S-phase start-time distribution. In other words, we must keep track of the laboratory clock compared to the cell-cycle clock.

Rate-equation approach
Our approach for calculating replication kinetics is based on a set of three coupled rate equations that describe the average space-time evolution of replication variables such as the replication fraction, f (x, t), or the replication fork density, ρ(x, t) (defined as the average number of forks per length at location x and time t). As we will see, the equations can be numerically solved in order to obtain the mean-field replication kinetics of a given replication scenario.
The first equation of our set describes the rate of change of the probability that a given position along the genome is replicated as a function of time, where the ± subscripts indicate the fork propagation direction (for right-and left-moving forks, respectively). Note that Table I presents the definition and units of all symbols used in this article. Equation 1 simply states that the rate of change of the replication fraction is given by the product of the local fork densities times the rate at which each fork replicates. Note that the velocity of leftand right-moving forks can be different in our formalism, although here they are kept equal. The second and third equations of our set express the space-time rate of change of the right-and left-moving fork densities as and respectively. The left-hand sides of these equations describe the change of fork densities due to their displacement along the genome. The right-hand sides present the contributions associated with initiations and coalescences of forks. The first term is simply given by the product of the local initiation rate, I(x, t), normalized by the probability that the genome is not already replicated at position x . The second term gives the frequency at which converging replication forks meet at x and t. This coalescence rate is proportional to the local densities of both types of forks and their relative speed when they merge. This rate must also be normalized by the probability that the local DNA is not already replicated, For clarity, we did not explicitly specify in Eqs. 1-3 that f (x, t), ρ ± (x, t), v ± (x, t) and I(x, t) can all change with both the position along the genome, x, and the time elapsed from the beginning of S phase, t. Given a replication scenario for the initiation rate, I(x, t), and the fork propagation speed, v ± (x, t), Eqs. 1-3 can be numerically solved for f (x, t) and ρ ± (x, t). The initial conditions must be set to f (x, 0) = 0 and ρ ± (x, 0) = 0 for all positions along the molecule. Such a numerical analysis can be carried out using various approaches. We developed and ran our own numerical code using IGOR Pro [36]. In the code, we discretize the modeled space-time region and approximate the partial derivatives in the equations describing the replication kinetics by finite differences. Knowing the initial conditions for all positions at time t = 0, one can then numerically integrate the system as a function of time [details about such numerical methods can be found in Ref. [37], for example]. The results presented in this article were obtained using the single-step explicit Euler method. We used dx = 0.002 Mb and dt = 0.01 hr (note that we must have dx/dt ≥ v in order to solve the system). For simulating or fitting to larger regions of the genome, the added efficiency of a more sophisticated numerical method for integrating the partial differential equations could be worth the extra effort to implement.

Boundaries and finite-size fragments
When modeling a whole genome using our calculation technique, one must specify the boundary conditions. The modeling we have described so far applies to linear DNA where replications forks vanish when they reach the ends of the molecule. Our simulation algorithm and calculation method can easily be adapted for the modeling of circular DNA molecules by using Periodic Boundary Conditions (PBCs). This simply means that any fork that reaches one end of our modeled domain (i.e., x = 0 or x = L) is re-inserted at the other end. In terms of our rate-equation system, using PBCs means imposing the conditions that ρ ± (0, t) = ρ ± (L, t) at all times t.
Molecular-combing experiments are often designed to extract DNA fragments from a particular region of the genome [e.g. [11,38]]. In order to model such fragments, we have two options: 1. We can model the whole chromosome and extract the information about the region that represents the experimentally observed fragment.
2. We can model only the fragment of interest, adding boundary terms to include replication by forks that initiate outside the fragment and then enter it during S phase.
The first option is a straightforward application of the model presented in Section and has been applied, for example, to the study of Chromosome VI in yeast (S. cerevisiae) [39]. However, if the chromosome is much longer than the fragment of interest, as is usually the case, it will be experimentally and computationally inefficient to study and model the whole molecule. In such a case, the second option can be an excellent approach, as long as we can model the effective flux of incoming forks at the modeled domain boundaries. The simplest situation is that no forks from outside the domain cross the boundaries. However, such a situation is realistic only when the "outside origins," if they exist, are far from the fragment's boundaries. More generally, we can include the contribution of origins near the fragment region in our model by replacing the periodic boundary conditions with time-dependent functions at either side that represent the flux of forks coming from outside the domain. The fork densities at the fragment boundaries are then [34] where The functions β ± (t) characterize the injection of right-and left-moving forks at x = 0 and x = L, respectively. For example, if the initiation rate outside the left and right edges is constant (I ± , respectively) in the semi-infinite region beyond the analyzed fragment, and if the fork velocity is constant as well, then β ± (t) = I ± t [34]. However, the real advantage in focusing on the fork fluxes ρ ± is that we do not need to have a detailed scenario of replication outside the region where data are collected but we only need to find a form for the time dependence of ρ ± that fits the data well. If it turns out that ρ ± = I ± t(1−f ± ) fits the data well, then we can assert that origin initiation outside the region, insofar as it affects the region under study, is equivalent to a zone of constant intitiation rate. In practice, we have found that the functional form ρ ± = I ± t(1−f ± ), with I ± as fit parameters, works well, because the closest origin or origin cluster to the boundary tends to dominate.

Post-treatment of the rate-equation solution
Once the solution of our rate-equation system is obtained for a given replication scenario [I(x, t), v(x, t), boundary conditions, ...], we can use it to calculate the space-time distributions of initiations and coalescences. These probability profiles tell us where and when such events are more likely to occur, deepening our knowledge about the replication kinetics.
The initiation density function is given by where is the average number of initiations per replication cycle between x = 0 and x = L. The initiation density function is thus simply a normalized version of the initiation term in Eqs. 2 and 3. Similarly, we can normalize the coalescence term to obtain the space-time probability density of observing a coalescence as with the average number of coalescences per cycle given by Other useful replication-related information can be derived from our rate-equation solution. For example, we can use f (x, t) and ρ ± (x, t) to derive the distributions of the first initiation or the last coalescence times [see Ref.

Model validation
We now compare our rate-equation solution to simulation results. The origin-initiation scenario we chose is presented in Fig. 2(a). It consists of a Gaussian-shaped initiation zone at the center of a 1 Mb fragment. We used the boundary conditions derived at the end of Section , i.e., a constant level of initiation outside the fragment region [parameters for I(x, t) are presented in caption of Fig. 2]. In order to simplify the interpretation of the results, we assumed that the initiation profile of Fig. 2(a) does not change with time, although our modeling formalism allows such variation. Similarly, we used a constant fork velocity (v = 0.15 Mb/hr = 2.5 kb/min). Again, our formalism also works equally well if the fork velocity depends on space and/or time. Figure 3 compares the calculated replication fraction and the initiation and the coalescence probability densities with simulation data. Part (I) of that figure helps to understand the physical meaning of our rate-equation solution. The figure illustrates that the calculated replication profile corresponds to the average profile we would get were we able to perform an infinite number of simulation iterations (the mean-field solution). Figure 3 shows that our calculated initiation and coalescence probability density functions agree with simulation results as well. More specifically, we used Eqs. 6 and 8 to calculate the average numbers of initiations and coalescences observed within the modeled domain per replication cycle (we obtained N i = 0.91 and N c = 0.84). These calculated values agree, within statistical fluctuations, with the simulation averages of N i = 0.91 ± 0.04 and N c = 0.82 ± 0.04, based on 500 simulated replication cycles.
Finally, note that although the number of simulated replication cycles is large, the average profiles presented in part (b) are still relatively noisy. Using rate equations gives us access to the exact average replication kinetics, whereas a simulation always gives a noisy approximation. Reducing the noise to usable levels can take a large number of simulations. On the other hand, simulations are good for understanding the level of statistical noise in The gray scales go from 0 (white) to 3.6/Mb/hr (black) for initiations and from 0 (white) to 1.8/Mb/hr (black) for coalescences. Both the simulation and the calculation were done using the initiation function I(x) presented in Fig. 2(a) (constant throughout time). The fork velocity was set to v(x, t) = 0.15 Mb/hr. a given experiment (by simulating the exact number of molecular fragments collected, etc.) and can also give information about more subtle features, such as correlations in the fluctuations of different measurements.

Modeling a synchronous combing experiment
As described in Section , our Monte Carlo algorithm can mimic a simple combing experiment of synchronously replicating cells. The gray bands in Fig. 4(a) present simulated combing results obtained using the replication scenario described in Fig. 2(a). These bands represent the average simulated combing profiles for various incorporation times of the nucleotide analogs (from 1 to 8 hrs). For each band, 500 combed molecules were simulated; the gray bands shown in Fig. 4(a) are the average replication profiles ± the standard error of the means.  Fig. 3(c-I)). (b) Similar data for asynchronous replication cycles. The start-time distribution was the same as in Fig. 5). The simulated data were collected as in part (a) except that the start times of the cycles were randomly chosen. The rate-equation solution was modified to account for the start-time asynchrony (see text and Fig. 5 for details).
These simulation results can be directly compared to the calculated mean-field replication fraction, f (x, t i ), where t i is the analogs supplementing time. The calculated combing profiles are presented as solid lines in Fig. 4(a) and agree with the simulated ones within statistical limits.

Correction for cell-cycle asynchrony
In the previous section, we demonstrated that our rateequation approach is consistent with averaged combing results when all cell cycles are perfectly synchronous. We now show how our calculated replication fraction f (x, t) can be modified in order to take into account the asynchrony of S phase start times in the population of measured cells. We consider two cases: (1) the distribution of start times is known, and (2) cell cycles are completely asynchronous. The distribution is Gaussian, with standard deviation σstart = 1.5 hr). Right: Transformed replication fraction profile, fasync(x, t) (see Eq. 9). As in Fig. 3, the dashed lines indicate the 1%, 40%, 70% and 90% contour curves.

S phase start-time distribution
Let t refer to the laboratory clock and let φ start (t) be the S phase start-time probability density function. If φ start (t) is known (or assumed), we can convolute φ start (t) with our rate-equation solution for f (x, t) to find the asynchronous replication curve: where f async (x, t) is the average observed replication function. The convolution operation is illustrated in Fig. 5 for our replication scenario solution example of Fig. 3(c-I) and a Gaussian distribution for the cell-cycle start times. The situation depicted here could correspond to an experiment where cells are imperfectly synchronized by drug treatment [see Ref [20] for an example]. In such a case, Eq. 9 can be used to correct our calculation.
The impact of an asynchronous start-time distribution on combing measurements is presented in Fig. 4(b). For this figure, the same simulation procedure as in part (a) was followed except that a Gaussian distribution with a standard deviation of 1.5 hr was used to adjust the cellcycle clock with respect to the laboratory time. Again, the gray bands represent the average and standard error of the mean obtained from 500 simulations, while the solid lines represent calculated profiles, f async (x, t i ).

Perfectly asynchronous cells
Equation 9 allows us to study a great number of partially synchronous cases but must be modified when the spread of start times becomes wider than the typical replication time. An extreme case is the limit of perfect asynchrony, where molecular fragments are collected at any time during the cell cycle, with a uniform probability distribution. In this limit, f asynch (x, t) is constant for all x and t and no useful information about the replication kinetics can be inferred. Norio et al. introduced another way to analyze data obtained for this case of perfect asynchrony [9]. The method uses only information from fragments that show evidence of active replication (i.e., that have at least one replication fork). The main difficulty is that the length of S phase-the time between the first initiation and the last coalescence-is different for each cell because initiations are stochastic and vary from cell to cell. Although the analysis is more involved and beyond the scope of this paper, it is still possible to extract replication parameters (initiation rates, fork velocities, etc.). See Ref. [34] for details.

Fitting data
An analytical formalism such as that presented in this article makes fast predictions of statistical profiles from a given replication scenario. Using this mathematical approach, we can infer replication parameters from experimental data more quickly and precisely than using simulation-based fitting methods.
To illustrate these possibilities, we performed a simple fit of the synchronous simulation data presented in Fig. 4(a) using the initiation profile as a fitting variable. More precisely, we fitted for the three parameters, {I 0 , µ, σ}, used to describe the Gaussian initiation peak of Fig. 2(a) and for the outside contribution, I − , as well (the left side contribution, I + , was set to zero since a free fit gave a value consistent with zero). The eight experimental curves of Fig. 4(a) were simultaneously fit using a standard routine [36]. The best-fit initiation profile is presented and compared to the original one in Fig. 2(b). Adding parameters, we can also infer the fork velocity, v, and the standard deviation of S phase start times, σ start (for asynchronous data). Note the need to account for correlations in the noise fluctuations to evaluate correctly the statistical errors of the fit parameters [see Section and Ref. [34] for details].
Fork density profiles DNA replication patterns, which depend on chromatin context, origin efficiencies and replication fork rates, are highly inhomogeneous both within the genome and across different tissue types. Knowledge of these patterns and how they evolve is crucial to understanding regulation of the cell cycle in both normal and diseased cells. As we discussed above, molecular-combing measurements provide replication fractions, f (x), at a given time t during S phase, and we have shown here that such profiles can be reproduced using our calculation technique. In addition to the replication state of the molecule, combing experiments also indicate the location of replication forks, and hence their density, at the time of nucleotide-analog incorporation. This additional information about fork densities can also be compared with our rate-equation solution. Thus, we could improve the inference technique presented in Section by simultaneously fitting the replication fractions and the two fork density functions using a common set of parameters. However, such "global fits" are not trivial. The correlations between fluctuations in different parts of the data-for example, how replication fractions correlate with fork densities-must be taken into account [34].

Fork arrest and DNA damage
Tumor cells are characterized by extensive karyotype alterations and aberrant DNA replication kinetics associated with deletions and amplifications. The altered kinetics has been found to depend on amplicon structure [33,40], and recently a Werner syndrome cell line was shown to exhibit altered rDNA patterns and elevated densities of unidirectional replication forks. The relationship between genetic structure and replication kinetics is poorly understood, but repeat-rich sequences in heterochromatin and elsewhere in the genome have an important impact on replication kinetics and genome stability [41]. Our rate-equation techniques can be easily adapted to study the impact of stochastic fork arrests on the replication kinetics. Using the model presented in this article, we can model the effect of a change in speed of the replication forks as they move along the genome, via the velocity profile v(x, t). Such an approach could be used to model zones where forks slow (for example, the replication slow zones discussed in [42]). The same approach can also be used to model a local or global change in fork velocity due to a checkpoint response. On the other hand, phenomena such as fork blocks due to DNA damage [43], where the fork velocity is different in each cell, need to be modeled differently. Our model can be adapted to such situations. For example, we have considered cases where forks stall at randomly located DNA lesions and may eventually restart after repair [34,44]. The main modification needed is to introduce two new fork densities, ρ ± (x, t), which represent the densities of left-and right-moving forks that have stalled. The rate of change of these new densities is governed by the stall rate of regular forks and the repair frequency of stalled forks. The replication kinetics would then be described using five differential equations instead of three. Two extra terms must also be added to the right-hand sides of Eqs. 2 and 3, to represent the stall and repair rates.

CONCLUSIONS
In this paper, we presented a simple and versatile simulation technique to model inhomogeneous DNA replication kinetics. We demonstrated that our modeling approach can be used to analyze or simulate molecular combing data as well as microarray-based results. We showed that our model can also be used to infer information about the replication kinetics from fitting typical experimental measurements.


Using Simple NLP Tools to Trace the Globalization of the Art World

We introduce a novel task, that of associating relative time with cities in text. We show that the task can be performed using NLP tools and techniques. The task is deployed on a large corpus of data to study a specific phenomenon, namely the temporal dimension of contemporary arts globalization over the first decade of the 21 st century.

Introduction
This paper reports on the use of fairly simple Natural Language Processing (NLP) technology as a tool for social research. We seek to understand the globalization of contemporary art, relying on releases for art events worldwide over the period of 1999-2012. A first phase of the project exploited Named-Entity Recognition (NER) to extract cities named in the corpus.
In this second phase of the research, we turn our attention to a novel task: the temporal dimension captured in these texts. By identifying the timing of geographically-cited events, we are able to examine how the history and projected future of the art world evolve alongside, and as a result of, its current geographical structure. To determine whether cities named in press releases refer to events that occur in the past, the present, or the future, we need methods for resolving time expressed in text relative to the time of the release of the text. We use the Stanford Temporal TaggerSUTime (Chang and Manning, 2012), as well as rules we have built on top of the Stanford part-of-speech tagger (Toutanova et al., 2003), to identify the temporal referent of each city mentioned. The two systems in combination perform slightly better than either does alone, and at a high enough level of accuracy to produce meaningful data for network analysis. We start by describing the project in more detail in Section 2, as well as the data we use in Section 3. Section 4 discusses our method for identifying the temporal location of named events. The networks we build using this data afford some preliminary insights into the dynamics of contemporary arts globalization, which we present in Section 5. Our aim here is not, however, to offer definitive conclusions about change in the art world. Rather, we present these analyses as proof of concept to demonstrate how a novel NLP task can help advance a particular agenda in social research.

The Goals of the Project
The shifting geography of contemporary art is of interest to social scientists as an instance of cultural globalization (Bourdieu, 1999;Crane, 2008). Scholars and contemporary art practitioners alike have remarked on changes in the geography of globally-significant art activity in recent decades (Vanderlinden and Filipovic, 2005;McAndrew, 2008;Quemin, 2006). Accounts of these dynamics often converge: for example, many have pointed to the role that speculation in the art of the developing world plays in fragmenting the once geographically-concentrated field of contemporary art (e.g., (Lee, 2012;Stallabrass, 2006)). Such claims remain largely a matter of conjecture, however, having never been subjected to systematic empirical scrutiny.
Because of the difficulties of compiling appropriate data, contemporary arts globalization has rarely been studied using quantitative methods. The few quantitative analyses that exist focus on particular cities (Velthius, 2013), individual artists (Buchholz, forthcoming), or a particular category of art-related events, for example fairs (Quemin, 2013). Our project relies on NLP techniques to solve the problems of data collection that have limited previous research endeavors. We extract information from a large and growing corpus of press releases for contemporary art events occurring worldwide. An initial phase of the project revealed for the first time an object that has resisted empirical description: the contemporary global art world. The research was able to track changes in the art worlds geographic structure over the period of 2001 until 2012. Named Entity Recognition (NER) was used to identify cities mentioned in the body text of announcements as the location of other significant contemporary art events. Pooling documents over a given time window yielded a dynamic network of cities, enabling the researchers to track the evolution of the art world. The analysis revealed that a distinct central core of cities all in Europe and North America enjoyed a remarkably stable ranking as top art locales over the period; however, as the art world expanded, connection to these capitals became less important as an entry to the art world for locations outside the global north (Rule and Brandt, 2013).
"The city that I believed was my past, is my future, my present; the years I have spent in Europe are an illusion, I always was (and will be) in Buenos Aires" (Borges, 1978). Borges' remarks capture a paradox about places of significance: important places are places where important things have happened and will happen-but this is so from the viewpoint of the present, which is always evolving. This insight implies that in order to understand the dynamics shaping the geography of significant art activity, analysts need to disentangle the real time in which attributions of significance to locations occur, from the relative timing of the events that make them significant. The current phase of our research attempts to assess how the geographical evolution of the art world shapes contemporary arts history and its projected future. It does so by considering how events mentioned in the corpus are located in time, relative to the event in question. We build upon the NER-based strategy developed in Phase 1 for identifying cities as the location of important events. Here, however, we distinguish cities as the location of events in the past, present, or future.

Data
We draw on corpus of press releases for contemporary art events worldwide, distributed over the premier email digest service for contemporary art. The digest serves a professional clientele, and sends out 3-5 press releases a day. These emails are free to the digest's 90,000 subscribers (compare to Artforum's circulation of 35,000). Press releases are both written and paid for by the institutions sponsoring the events they announce. The corpus covers the years 1999 to 2012; it contains 10,566 documents and 6362284. Our automated detection, checked by human annotators, identified 1007 unique cities.
We conceive of each press release as an attempt to raise the status of the event in question by mentioning other high-status events to which it is related. Pooling the references of cities across documents thus gives rise to a network, analogous to a citation network. The press releases are rather conventionalized, and name cities almost exclusively as the location of important related events. The format is not just used by this particular digest service, but is common to the art world generally. A sample is shown in Figure 1.

Time Identification
Cities associated with art events are identified using an NER-based approach, which performs at an f-measure of 64.5 and which we do not further describe in this paper. We use two approaches for temporal resolution to assess whether the cities mentioned are the locations of events that happened in the past, are currently happening, or expected to happen in the future.
The first approach analyzes explicit time ex-  Table 1: Accuracy and precision, recall and F-measure for past, current and future events. * Baseline: Tagging all events as current.
pressions in dates, time durations, holidays, etc. The second approach uses verbal tense to resolve grammaticized reference to time. In both approaches: (1) We use distance within a sentence as a heuristic to associate a temporal expression to a city.
(2) Cities associated with art events are tagged as past, current or future, where current is the default tag. This section describes both approaches and how they can be combined to enhance the performance.

Explicit Lexical Temporal Expressions
Explicit temporal expressions could be partial or complete dates (Jan-01, March 2nd, 2014), specific named days (Christmas 2001), deictic time expressions (last Wednesday), time durations (the past five years), seasons (next fall), etc. We use the Stanford Temporal Tagger, SUTime (Chang and Manning, 2012) to detect and evaluate Temporal expressions relative to the publication date of the press release. Temporal expressions that evaluate to incomplete dates take on the missing date components from the publication date. If the temporal expression resolves to a date before, same as or after the publication date, then the event associated with the city is tagged as past, current, or future event respectively. Cities that are not associated with explicit temporal expressions are given the default tag: current.

Grammaticized Temporal Expressions
Tense is grammaticalization of location in time (Comrie, 1985). Our tense-based approach has two steps. First, the Stanford part-of-speech tagger (Toutanova et al., 2003) is used to identify the POS tags of all tokens in the documents. Second, we use hand-written rules to identify the temporal interpretation of every event. The rules are mostly based on the POS tags of verbs; we use the tag sets from the Penn Treebank Project (Marcus et al., 1993). We use only the verbal tags: VB (Verb, base form), VBD (Verb, past tense), VBG (Verb, gerund or present participle), etc. and MD (Modal). Here are some examples of the rules; events associated with: (1) VBP or VBZ are tagged as current.
(3) VB that are preceded by will/MD are tagged as future.

Results and Combining of the Two Approaches
During development, we found that the two time tagging approaches perform well in different situations. For example, the lexical approach has the higher recall (66.1%) for past event as writers often mention the exact date of past events. However, in the case of current events dates are seldom mentioned, the present tense is enough to communicate that the event mentioned in a sentence is ongoing. This accounts for the higher recall of the grammaticized (89%) as against the lexical approach (75.5%) for current events. If we do not assign current by default to the failed cases in both the lexical and grammaticized approaches, the recall for current events drops drastically for the lexical approach (75.5% 0.69%) compared to the grammaticized approaches (89% 81.5%).
Combining the two approaches improves the performance of the time tagging task. We start with one approach to tag an event and when that approach fails (does not provide an answer), the other approach kicks in. If both approaches fail, the event is tagged by default as current. For example, in approach Grammaticized Lexical (G L ), when the Grammaticized approach fails to tag an event, the Lexical approach is used. The best combination is achieved by running the Lexical Grammaticized (L G ) first; if the output tag is past, accept its answer; otherwise run the Grammaticized Lexical system (G L ) and accept its verdict. If G L has no answer, choose present. The accuracy of this approach is 78.0%. For more details, see table 1. Please note that analyses in section 5 rely on the Grammaticized Lexical (G L ) approach.  Table 2: Properties of networks of past, present, and future art events, for three periods. "Share" refers to the percentage of events in relative time during a given period.

Creation and Interpretation of Networks
To analyze the geography of the art world, we generated networks of events in relative time over three four-year periods: 2001-2004; 2005-2008; and 2009-2012. Nodes represent cities named as the location of art events, and edges co-citations from a particular geographical and temporal standpoint. Specifically, we attribute ties between cities of the past, present or future of a given period when they are mentioned together in press releases for events currently occurring in the same city.
Ties are weighted to reflect the number of comentions. For example, if Basel is noted twice as a future event location by press releases for events in New York, and LA is mentioned by New York press releases three times, the edge between Basel and LA receives an edge weight of 2, reflecting the number of shared mentions. Basic properties of the nine resulting networks are displayed in table 2. Notice the marked difference in the distribution of events in relative time between 2005-2008 and 2009-2012. Cities mentioned in connection with events in both the future and the past figure as a greater proportion of the total in the last, as compared to the middle period. Though we consider the possibility that this is an artifact of the data, it seems more to likely that the shift reflects the impact of a real world development: namely, the global economic recession, which began in 2008 and had profound effect on funding for contemporary art. In the context of the economic contraction, it seems that both contemporary arts future and its history become more important relative to its present.  Recall that the graph captures the world's future, as it is structured by the geography of significant art world activity during this period. Just as in Borges remarks about Buenos Aires, we observe here how important places in the art worlds present shape the map of its future.
This graph suggests an art world that is rather conservative, more so, at least, than accounts of the speculation driving its globalization would imply. The network exhibits a classic core-periphery structure; the cities at its center are the "big" art capitals of the global north: London, New York, Berlin, Paris, Venice. In other words, the common proposition that the "hype" around emerging art markets has de-centered the contemporary art world is not borne out by a close empirical examination from the perspective of the hype-generators themselves. Rather, it would seems that the role of these cities as the prime location for contemporary art in the present enables them to project themselves as the location of significant art activity in the future.


extfieldsofstudy:
Computer Science
Medicine

text:
Selection of entropy-measure parameters for knowledge discovery in heart rate variability data

Background Heart rate variability is the variation of the time interval between consecutive heartbeats. Entropy is a commonly used tool to describe the regularity of data sets. Entropy functions are defined using multiple parameters, the selection of which is controversial and depends on the intended purpose. This study describes the results of tests conducted to support parameter selection, towards the goal of enabling further biomarker discovery. Methods This study deals with approximate, sample, fuzzy, and fuzzy measure entropies. All data were obtained from PhysioNet, a free-access, on-line archive of physiological signals, and represent various medical conditions. Five tests were defined and conducted to examine the influence of: varying the threshold value r (as multiples of the sample standard deviation σ, or the entropy-maximizing rChon), the data length N, the weighting factors n for fuzzy and fuzzy measure entropies, and the thresholds rF and rL for fuzzy measure entropy. The results were tested for normality using Lilliefors' composite goodness-of-fit test. Consequently, the p-value was calculated with either a two sample t-test or a Wilcoxon rank sum test. Results The first test shows a cross-over of entropy values with regard to a change of r. Thus, a clear statement that a higher entropy corresponds to a high irregularity is not possible, but is rather an indicator of differences in regularity. N should be at least 200 data points for r = 0.2 σ and should even exceed a length of 1000 for r = rChon. The results for the weighting parameters n for the fuzzy membership function show different behavior when coupled with different r values, therefore the weighting parameters have been chosen independently for the different threshold values. The tests concerning rF and rL showed that there is no optimal choice, but r = rF = rL is reasonable with r = rChon or r = 0.2σ. Conclusions Some of the tests showed a dependency of the test significance on the data at hand. Nevertheless, as the medical conditions are unknown beforehand, compromises had to be made. Optimal parameter combinations are suggested for the methods considered. Yet, due to the high number of potential parameter combinations, further investigations of entropy for heart rate variability data will be necessary.

Background
Heart rate variability (HRV) is the variation of the time interval between consecutive heartbeats. It highly depends on the extrinsic regulation of the heart rate (HR) and reflects the balance between the sympathetic and the parasympathetic nervous system [1]. Batchinsky et al. [2] develop a collection of methods for using HRV to describe regular periodic oscillations in the heart rate, attributed to the vagal and/or sympathetic branches of the autonomic nervous system.
Research on HRV has attracted considerable attention in the fields of psychology and behavioral medicine. It has its origin in the search for non-invasive correlates of injury severity which can be extracted from available signals in order to discover new cardiac biomarkers [2]. These signals are usually ones that are routinely measured, and include sources like a photoplethysmogram or the electrocardiogram (ECG) [3]. Electrocardiography is an interpretation of the electrical activity of the heart over some period of time. It is a non-invasive procedure, using electrodes attached to the surface of the skin, and is commonly used to measure the heart rate, the regularity of the beats, and characterize properties or injuries in the heart chamber. An R-to-R interval (RRI) describes the latency between two consecutive R peaks in the ECG. The RRI time series are used as input for the determination of HRV parameters.
In studies of HRV, both time-and frequency-domain measures are typically used by practitioners and researchers [1,4]. Additionally, further knowledge about the subject's status can be discovered by the evaluation of certain patterns and shifts in an "apparent ensemble amount of randomness" of a stochastic process [5]. This randomness, as well as the predictability of this process, can be measured by entropy [6]. Thus, it is a commonly used tool to describe the regularity of large biomedical data sets. The more regulatory inputs a system has, the higher its irregularity is, due to interference of those regulatory systems. This assumption is also true for many biomedical systems such as HRV [7]. Therefore, it is a reasonable hypothesis that a more regular heart rate variability is connected to a defect in a regulatory system. To measure this irregularity, Pincus et al. proposed approximate entropy (ApEn) [8]. Further types of entropy were developed based on this method to improve it [9]. Because of its ability to measure regularity, entropy is widely used as a diagnostic tool in medicine to derive and discover biomarkers in large biomedical data. Its applications range from sudden infant death syndrome [7], to complexity analysis of intracranial pressure dynamics during periods of severe intracranial hypertension [10], to quantification of amplitude variations in mechanomyographic signals [11], to analysis of short gait data sets [12], to automatic detection of normal, pre-ictal, and ictal conditions from recorded electroencephalography signals [13], to the postural sway in stroke patients [14].
The approach to quantify the structural complexity (or, inversely, regularity) of the HRV is called heart rate complexity (HRC) and utilizes methods derived from nonlinear dynamics. Note that complexity and variability are not necessarily the same [15]. A periodical signal, such as a sine wave, is variable but not complex. This property allows complexity measures to ignore the complicated periodic oscillations to at least some extent [16].
To date, numerous entropy types are widely accepted as measures of the HRC. However, since their function is controlled by three, four or six parameters, there are many possible combinations to choose from. The selection of criteria for these parameters is controversial and heavily depends on the intended purpose and the data at hand [17]. The variation of only one parameter often results in highly non-linear behaviour [12]. Therefore, results of calculations with different parameters cannot be easily extrapolated from existing data, but have to be computed individually. Hence, the variation of several parameters in order to optimize the reliability of the results is a very time consuming process. To be used in daily routine, reasonable parameters have to be selected prior to the evaluation of HRC.
The determination of appropriate parameters for entropy in general [11,12,[18][19][20] and for HRV applications in particular [16,17,21,22] is the subject of ongoing research. A summary of the current knowledge regarding parameter selection is given in the subsection "Parameter selection". In order to verify the results of previous publications and to extend them by the usage of more entropy measures, this work focuses on the parameter selection for approximate, sample, fuzzy and fuzzy measure entropy and their implications for HRV data. The objective of this paper is to provide a reference for choosing parameters for HRV applications based on their influence on the different entropies' ability to distinguish significantly between pathological and non-pathological recordings. The main research questions addressed in this paper are: (1) what does the choice of the threshold value(s) mean for the entropies to be a direct measure of regularity, (2) how does the data length influence significance for different data sets, (3) how should the weighting factor(s) be chosen for fuzzy and fuzzy measure entropy, and finally, (4) what are the constraints for choosing the threshold values?.

Methods
In the following sections approximate (ApEn), sample (SampEn), fuzzy (FuzzyEn) and fuzzy measure entropy (FuzzyMEn) are described in detail. They are built on each other and ordered by increasing complexity. The number of parameters increases from three for ApEn to six for FuzzyMEn, with the basic parameters staying the same and being extended by new ones in each step. Afterwards, the challenge of parameter selection and the tests conducted for the parameters to the different entropy types are described. Finally, the data used for all performed tests are briefly described.

Approximate entropy
The main idea behind approximate entropy is that a sequence is regular if a subsequence and an expansion of the subsequence are similar. It was developed by Pincus [8] and is calculated the following way.
Given a sample sequence {u 1 , . . . , u N }, a template length m and a threshold value r, the sequence is first split into overlapping sequences Next, define C m i as the number of j = 1, . . . , N − m + 1, is the Chebyshev distance, i.e., the maximum distance between two elements of X m i and X m j . Constructing the same sequences, but with template length m + 1, yields C m+1 i . Then, φ m and φ m+1 are defined as: The approximate entropy is then defined as

Sample entropy
Richman and Moorman showed in [23] that approximate entropy is biased towards regularity. Thus, they modified it to sample entropy. The main difference between the two is that sample entropy does not count self-matches, and only compares the first N − m subsequences instead of all N − m + 1, so the same amount of subsequences are used in j m and j m+1 [23]. It is calculated in the following way.
Given a sample sequence {u 1 , . . . , u N }, a template length m, and a threshold value r, first the overlapping sequences, {X m 1 , . . . , X m N−m } are constructed as for ApEn. As opposed to ApEn, C m i is now defined as the number of j = 1, . . . , is again the Chebyshev distance. Applying the same for template length m + 1 results in:

Fuzzy entropy
ApEn and SampEn are very sensitive with respect to the threshold parameter r. They show, on one hand, a very abrupt behavior, and on the other hand less significance for small r, as was discussed in [9] and can be seen in figure 1. To soften these effects, Chen et al. developed fuzzy entropy, which uses a fuzzy membership function instead of the Heaviside function [9]. FuzzyEn is calculated as follows.
Let m and r again be the template length and the threshold value and n a weight for the fuzzy membership function. Sequences {X m 1 , . . . , X m N−m+1 } are defined as for ApEn and SampEn from an input sequence {u 1 Next, using a fuzzy membership function x µ(x, n, r) , n, r) . According to Chen et al. [9], fuzzy membership functions must be continuous and convex. The first property guarantees only slow similarity changes, and by the second condition self-similarity is a maximum of the function. For all tests conducted in this study, x → e −(x/r) n was chosen as the fuzzy membership function as in [9]. For n ∞ this fuzzy membership function converges to the Heaviside function. These steps are repeated for template length m + 1 to get D m+1 . Consequently, φ m and φ m+1 are calculated the following way: Fuzzy entropy is defined as FuzzyEn(m, r, n) := lim and can be estimated by FuzzyEn(m, r, n, N)

Fuzzy measure entropy
Liu et al. proposed adding a function for global similarity to the fuzzy entropy and called the combination fuzzy measure entropy [24]. It can be calculated as follows. Given a data sequence {u 1 , . . . , u N }, a template length m, two threshold values r L and r F , and two weighting parameters n L and n F (r L and n L correspond to the local term and r F and n F to the global one), a sequence {X m 1 , . . . , X m N−m+1 } is constructed like before. In the next step, it is transformed into a local sequence (5) and XF m i := X m i − u mean , and u mean is the mean value of the complete data sequence {u 1 , . . . , u N }.
Using the local and the global parameters for the fuzzy membership functions, the matrices DL m and DF M are defined as DL m Afterwards, all these steps are repeated for template length m + 1 to get DL m+1 and are defined as: Fuzzy measure entropy is then defined as FuzzyMEn (m, r L , r F , n L , n F ) := lim

Parameter selection
As one can see in the description of each entropy, the various entropy types have three, four or six parameters. Thus, there are many possible combinations to choose from. The parameters which are varied in the test cases, their ranges and values mentioned in the literature, and the choice of certain fixed parameters are described here. The most common choice for the template length is m = 2, as it was recommended by Pincus and Goldberger for ApEn [7], by Yentes et. al for SampEn [12], and confirmed by other studies, e.g. [20]. A false nearest neighbor method is also sometimes used, but according to Chon et al. the standard choice leads to the statistically best solutions for ApEn for human heart rate variability data [25]. In [9,24], the template length is set to m = 2 for other entropy types as well. For comparability, this value was used for all tests.
Some publications describe a sensitivity of the entropies to the data length N [1, 6,8]. Therefore, the significance of the entropies calculated of parts of the heart rate variability data has been tested with increasing data set size.
For the threshold parameter r, Pincus suggested in [8] to choose a value between 0.1 s and 0.25 s, where s is the sample standard deviation of the data sequence. This is also the standard range in most publications, with the most common choice of r = 0.2 s [17,21,23,24,26]. To examine threshold parameters in the standard range, they were tested with heart rate variability data.
In [19,27], the so called flip-flop effect is described, where for some values of the threshold value r a signal has a higher entropy compared to another and for other choices of r a lower one. This is also shown in [21] for heart rate variability data. To test for this effect, various entropies of signals from one database were calculated with parameters inside the standard range. Lu et al. also showed this effect in [26]. Therefore, they proposed choosing r ∈ [0.1 · s, 1.0 · s] for ApEn in such a way as to maximize the entropy value. Since finding a maximum is computationally very expensive, Chon et al. created in [25] an empirical formula to calculate an r, hereinafter called r Chon , which approximates a maximizing r. For m = 2, it can be formulated as: where s 1 is the standard deviation of the distances in the data sequence, i.e., the standard deviation of {(u 1 − u 2 ), . . . , (u N−1 − u N )}, and s 2 the standard deviation of the complete data sequence. This formula was derived from nonphysiological data and Liu et al. showed in [17] that it is not always a good approximation of the maximizing r, but actually leads to more significant results than the maximizing r, when applied to heart rate variability data.
Regarding fuzzy function parameters, Chen et al. [9] used n = 2 and r = 0.2 s for test signals. They described in [9] that for a larger n, the closer data points are weighted more strongly. Liu et al. [24] used the weighting factors n L = 3, n F = 2 and r L = r F = 0.2s for heart rate variability analysis. Their choices for n L and n F were given without any motivation.

Test cases
To get further knowledge of the parameters, heart rate variability data were used to compare different choices of r L , r F and n, n L and n F with respect to the resulting significance of the statistical tests comparing pathological and normal cases.
The following test cases using the data described in the Data section have been conducted to answer the research questions stated in the Background section and to support an "optimal" parameter selection for all entropy types for heart rate variability data. For each research question, one test case has been designed (except for the third, which is covered in two test cases: one for each entropy type under consideration). All tests have been conducted consecutively. Fixed parameters have been taken from literature or based on the outcomes of preceding tests.

Test case 1
Variation of the threshold values r = r L = r F within the standard interval [0.1 s, 0.25 s] [8] to show its influence on the entropy values and the aforementioned flip-flop effect [21,27]. Fixed parameters were m = 2, n = n L = 3 and n F = 2 [7,24]. For data length N, the data length of the shortest RR interval sequence of the available data N = 1126 sets was used.

Test case 2
Variation of the data length N = x · 110 with x ∈ [1, 10] to show its influence on the significance. Maximum N was 1100 due to the length of the available test data. Fixed parameters were m = 2, n = n L = 3, n F = 2 and r = r Chon or r = 0.2 s [7,17,21,23-26].

Test case 3
Variation of the weighting factor n for FuzzyEn in the interval [1,6] to show its influence on the significance. Fixed parameters were m = 2 and r = r Chon or r = 0.2 s [7,17,21,[23][24][25][26]. N = 1000 was chosen based on the results of test case 2.

Test case 4
Variation of the weighting factors n L and n F for Fuzzy-MEn in the interval [1,6] to show their influence on the significance. Fixed parameters were m = 2 and r L = r F = r Chon or r L = r F = 0.2 s [7,17,21,[23][24][25][26]. N = 1000 was chosen based on the results of test case 2.

Test case 5
Variation of the threshold values r L and r F for Fuzzy-MEn in an interval of [0.25 · r Chon , 6 · r Chon ] and in an interval of [0.1 s, 0.25 s] [8] to show their influence on the significance. The parameter m = 2 was fixed [7]. N = 1000, n L and n F were chosen based on the results of the previous tests.
In order to test for their statistical significance, the calculated entropies were first tested for normality using Lilliefors' composite goodness-of-fit test [28]. If this test was positive for all results within a subset of a test case (i.e., certain database and/or choice for r), the p-value was calculated with a two sample t-test, and otherwise a Wilcoxon rank sum test was performed. To ensure comparability, the same statistical test was used for each subset of a test case.
Since statistical tests used in this work are based on the same null hypothesis, the same subject groups, the same endpoints and only slight variations of the analysis method, interaction of the observed results is not only possible, but highly probable. On the other hand, p-value adjustments such as the commonly used Bonferroni correction assume uncorrelated endpoints and are therefore considered inappropriate for the tasks in this work [29]. Besides, the aim of this work is not to test whether there is a difference between groups, but to investigate the ability of varied methods to detect those differences.
Due to the high computational complexity of the tests conducted (e.g., the calculation of the variation of the weighting factors n L and n F for FuzzyMEn takes several hours for only one r and one comparison of databases), the authors had to refrain from more robust randomized testing strategies (i.e., the permutation of the original data), since the computation time would multiply by at least a thousand times.

Data
All data used for the described tests have been taken from http://Physionet.org [30], a free-access, on-line archive of physiological signals. They are described in detail in this section.
To create a control group all databases described as non-pathological were combined into one database, afterwards called np. This database contains the Normal Sinus Rhythm RR Interval Database, which consists of the beat annotations of 54 long-term ECG recordings, digitized at 128 samples per second, of subjects with normal sinus rhythm (30 men, aged 28.5 to 76, and 24 women, aged 58 to 73). Furthermore, it includes the MIT-BIH Normal Sinus Rhythm Database, which consists of 18 long-term recordings (5 men, aged 26 to 45, and 13 women, aged 20 to 50) digitized at 128 samples per second. The Fantasia Database of 120-minute recordings of twenty young (10 men and 10 women; 21 -34 years old) and twenty elderly (10 men and 10 women; 68 -85 years old) healthy subjects with ECG digitized at 250 Hz was also used [31]. The record "fantasia/f1o09" had to be excluded due to its high number of supraventricular premature beats. This results in a total database size of 111 recordings. RR intervals greater than 2.5 seconds were excluded to ignore artifacts.
Two databases were used to search for pathological effects. One was the Congestive Heart Failure RR Interval Database, afterwards referred to as chf2, which includes 29 long-term ECG recordings, with a sampling frequency of 128 Hz, of subjects aged 34 to 79 (8 men and 2 women; gender not known for the remaining subjects) with congestive heart failure (NYHA classes I, II, and III) [32].
The second one was the MIT-BIH Arrhythmia Database, afterwards called mit, which contains 48 half-hour recordings, sampled with a frequency of 360 Hz, from 47 subjects (25 men aged 32 to 89 years and 22 women aged 23 to 89 years) [33]. It contains a set of randomly chosen signals and 25 signals especially chosen to include examples of uncommon but clinically important arrhythmias recorded at the BIH Arrhythmia Laboratory [33].
The two databases chf2 and mit were always evaluated separately to keep the results as homogeneous as possible and to avoid the mutual neutralization of ab-normalities. To ensure comparability, the data length N was equal for all recordings. Longer recordings were cropped at the beginning and the end in equal shares.
In a way, this work can be considered as a pilot study, as only previously recorded data are used. Furthermore, the findings of this work will be incorporated in followup studies.

Results
The following sections show the results of our tests, which are described in the Test cases section above. An overview of these results is given in Table 1.
The results were tested for normality using Lilliefors' composite goodness-of-fit test. Throughout the whole test case 2, which combines results of two different ways to determine r (see figure 2), data were either normally or not normally distributed. Hence, entropy values were compared using the Wilcoxon rank sum test. Entropies calculated with r as multiple of s were distributed normally in our test scenarios and could therefore be compared using the two sample t-test (as in figures 3 (B, D), 4 (B, D) and 5 (B, D)). Entropy values derived using r as multiple of r Chon , on the other hand, were not normally distributed and hence compared using Wilcoxon's rank sum test (as in figures 3 (A, C), 4 (A, C) and 5 (A, C)). Only findings of the same statistical test are combined in the following sections. Figure 1 corresponds to test case 1 and shows the effect of different threshold values on the entropy values. One of the defining characteristics of FuzzyEn and FuzzyMEn, their relative insensitivity to changes in r, is clearly visible. The flip-flop effect can be observed for ApEn. For the chf2 database, entropy values are higher for pathological data where r <0.15 s, whereas they are lower for r >0.15 s. This behavior is reversed for the mit database, with lower entropy values for pathological data for r <0.18 s and higher entropy values for pathological data for r >0.18 s.
The results for test case 2, analyzing the sensitivity of the entropy types to the data length are represented in figure 2. The more data points evaluated, the higher the Variation of the weighting factors n L and n F for FuzzyMEn in the interval [1,6] Best results with n L = 2, n F = 1 Best results with n L = 1, n F = 3 separation between pathological and non-pathological data. When using the threshold value r = 0.2 s, significance is already reached with N ≥ 200 data points, whereas with r = r Chon more data (N ≥ 1000) are needed before significance is reached. Comparing the different methods when using r = 0.2 s, one can see that they only differ when N is very small. With higher N , their behavior converges. Figure 3 presents the results of test case 3, showing that the results become insignificant for n >3 with r = r Chon , and for n >1 with r = 0.2 s in the mit database (C, D). The increasing p-values are below the significance level (p <0.05) for the chf2 database within the observed range of n (A, B). , calculated with n L = 1 and n F = 3. In (A), the best results are achieved with r L < ∼ 1 · r Chon or r F < ∼ 1 · r Chon , whereas they have to be greater than or equal to 1 · r Chon in (C) for good performance. In (B), r L <0.12 s or r L > ∼ 0.2σ yield significant results, however rF does not matter that much. In contrast, the best performance in (D) is achieved with r F > ∼ 0.18σ and r L ∈ [0.14 s, 0.22 s].

Discussion
As a summary of the following discussion of the previously presented results, the combinations of parameters which yielded the best results are listed in Table 2. Some of the tests concerning the parameter selection showed no clear results. A couple of evaluations revealed contradicting results depending on the database (chf2 vs. mit ) or on the chosen parameters (e.g., r = r L = r F = r Chon vs. r = r L = r F = 0.2s). The latter is easy to overcome by choosing a different set of parameters based on the choice of the threshold value r. However, a compromise must be found in order to find parameters suitable for both databases.
One of the hardest difficulties lies in the choice of the threshold value r due to the flip-flop effect, i.e., for some parameters one data set has a higher entropy compared to another, but this order is reversed for different parameter choices [17,27]. This can occur for simple signals, but also when analyzing heart rate variability data, as in figure 1. This leads to difficulties with the interpretation of the entropy, i.e., the direct assignment of entropy values to pathological or non-pathological data without a given r. In our tests, the effect occurred for ApEn, but it is also reported to occur for SampEn and FuzzyEn as well by Boskovic et al. [27]. The flip-flop effect does not allow us to make a clear statement, as in [7][8][9]17], that a higher entropy corresponds to a higher irregularity. Since this effect can happen for all different entropy values [27], they should not be seen as a direct measure of regularity, but rather as an indicator of differences in regularity with regard to certain time periods.
Compared to the threshold value, the data length N seems to have a smaller effect on the ability of the entropy measures to differentiate pathological from non-pathological data sets. Generally, a larger N leads to a higher probability of significance. If possible, the data length N should be longer than 200 data points when using r = 0.2 s. This finding is consistent with previous studies, e.g. [12]. Surprisingly, it should even exceed a length of 1000 data points when using r = r Chon , assuming a continuing trend in figure 2. Due to the length of the available test data, the range had to be restricted for this test. In case of HRV data, N is the number of recorded heartbeats and therefore proportional to the duration of the recording. Thus, to increase N, longer measurements are necessary. The Task Force of The European Society of Cardiology and The North American Society of Pacing and Electrophysiology [4] recommends a duration of five minutes for short time recordings, which would result in 300 data points at an average heart rate of 60 beats per second. This would be sufficient when using r = 0.2 s, but not for r = r Chon . These considerations should be kept in mind when dealing with HRV data.
Due to the different behavior when varying n, n L and n F given different threshold parameters r = r L = r F = r Chon and r = r L = r F = 0.2 s, the parameters n, n F and n L have been chosen independently for the different threshold values. This is no constraint to the method, since the choice of r is known beforehand anyway (and is not based on the potentially unknown medical condition of the subject). The values n L = 2 and n F = 1 for r L = r F = r Chon , and n L = 1 and n F = 3 for r L = r F = 0.2 s showed better results than n L = 3, n F = 2 as proposed by Liu et al. in [24]. Unsurprisingly, similar values were found for n, as n for FuzzyEn equals n L for FuzzyMEn. For consistency, n = n L = 2 for r = r Chon and n = n L = 1 for r = 0.2 s are recommended.
The tests concerning r F and r L showed that there is no optimal choice, since the results in figure 5 (A) and (C) as well as (B) and (D) contradict each other. Nevertheless, both r F = r L = r Chon and r F = r L = 0.2 s, as described in the literature [9,24,25], are the most reasonable compromise between figure 5 (A) and (C), and (B) and (D), respectively.
Finally, a number of important limitations of this study need to be considered. First, this study is limited to previously recorded signals of only two different cardiac diseases due to the availability of data. As already mentioned in the Data section, a separate evaluation is warranted, to avoid the mutual neutralization of abnormalities. Furthermore, the template length was fixed to m = 2 for all calculations, as the number of possible variations of parameters would get too high otherwise and m = 2 seems to be a reasonable choice in all investigated literature, e.g., [8]. In [34] it was reported that spikes due to recording errors in heart rate variability data can disturb ApEn and SampEn. No tests were done to examine this behavior for the included entropies. Finally, the study did not evaluate the dependency of the entropies on age and gender as reported in literature [35,36].

Conclusions
The results of the presented study clearly stress the need for further investigations of signal entropy for heart rate variability data. Given the wide range of different medical conditions of subjects, the assortment of available methods to calculate entropy, and their customizability with up to six degrees of freedom, it is almost impossible to cover all combinations in a single study. Future work will therefore be focused on overcoming the limitations of the presented work, i.e., extending the evaluations to other cardiac diseases, the variation of the template length m, investigating the robustness with respect to recording errors, and the relation of the parameter choice to age and gender.