Ali Makke

Osman Salem

Mohamad Assaad

Yvette Hassine Moungla

Ahmed Mehaoua

Flooding Attacks Detection in Backbone Traffic Using Power Divergence

Keywords: C.2.3 [Network Operation]: Network Monitoring Algorithms, Design, Experimentation, Measurement, Performance, Security Network Anomaly Detection, DDoS, Power Divergence, Hellinger Distance, Sketch

Flooding attacks detection in traffic of backbone networks requires generally the analysis of a huge amount of data with high accuracy and low complexity. In this paper, we propose a new scheme to detect flooding attacks in high speed networks. The proposed mechanism is based on the application of Power Divergence measures over Sketch data structure. Sketch is used for random aggregation of traffic, and Power Divergence is applied to detect deviations between current and established probability distributions of network traffic. We focus on tuning the parameter of Power Divergence to optimize the performance. We evaluate our approach using real Internet traffic traces, obtained from MAWI trans-Pacific wide transit link between USA and Japan. Our results show that the proposed approach outperforms existing solutions in terms of detection accuracy and false alarm ratio.

INTRODUCTION

Security threats for computer network have increased significantly, which include viruses, exploitation of software vulnerabilities, worm-based attacks, Distributed Denial of Service (DDoS), etc. DDoS through TCP SYN flooding is able to make silent any web site, especially with the use of BOTNETs (roBOT NETworks) containing large number of slave machines (zombies).

Many web sites suffered from SYN flooding attack that aims to exhaust server resources and to deny access for legitimate users. Recently, PayPal has been driven offline after dropping WiKiLeaks donations (operation PayBack in 2010). The list of victim web servers is very long (Twitter, CNN, Yahoo, Amazon, Ebay, DoubleClick, etc.). Therefore, flooding attacks need to be detected in order to cope with ongoing anomaly as soon as possible.

With the distributed nature of DDoS attacks, and its impact on the performance of the routers, the detection and reaction mechanisms must be pushed to the core network (Backbone), or near the sources of attack. However, with the growing complexity in analyzing huge amount of data traffic in backbone network, the analysis of each traffic flows is unscalable and computationally expensive.

Many change point detection algorithms have been proposed, and applied to the time series resulted from the aggregation of whole network traffic in one flow. These methods are based on the identification of the change point where heavy deviation occurs in the resulted time series. However, as these methods aggregate the whole traffic in one time series, low intensity attacks is buried by the variation of background traffic in high speed networks, and may pass undetected. Moreover, the time series derived from IP traffic are subject to many variations that are irrelevant to anomaly. The time series are non-stationary and tend to change over time, leading to a lot of false alarms. Reducing the false alarms and increasing the detection accuracy in such methods are a challenging problem.

In this paper, we propose a new framework for the detection of flooding attacks. The proposed approach is intended to detect low intensity attacks in backbone traffic. It is based on Power Divergence (PD [START_REF] Havrda | Quantification method of classification processes: The concept of structural α-entropy[END_REF]) and Sketch data structure [START_REF] Krishnamurthy | Sketch-based Change Detection: Methods, Evaluation, and Applications[END_REF].

Sketch is used for dimensionality reduction and to derive time series for randomly aggregated traffic, and PD measures the difference between 2 set of probability values, and detects the deviation between these sets in normal and under anomaly conditions. The PD is a general form of divergence measures, where Kullback-Leibler (KL [START_REF] Cover | Elements of information theory[END_REF]), Hellinger Distance (HD [START_REF] Tang | Sketch-Based SIP Flooding Detection Using Hellinger Distance[END_REF]), Chi-square (χ 2 [START_REF] Taneja | Bounds on Triangular Discrimination, Harmonic Mean and Symmetric Chi-Square Divergences[END_REF]), etc. can be derived from PD by changing its parameter (β). We emphasize in our analysis the optimal parameter of Power Divergence.

The rest of this paper is organized as follows. Related work is provided in section 2. Section 3 describes briefly the Sketch data structure. The proposed approach is explained in section 4. In section 5, we present our experimental results. Finally, in section 6 we present the conclusion and our future work.

RELATED WORK

Several approaches have been proposed for network anomaly detection, and they are based on different techniques, such as Haar-wavelet analysis [START_REF] Lu | Network Anomaly Detection Based on Wavelet Analysis[END_REF], entropy based method [START_REF] Lakhina | Mining Anomalies using Traffic Feature Distributions[END_REF], CUmulative SUM (CUSUM) algorithm [START_REF] Siris | Application of Anomaly Detection Algorithms for Detecting SYN Flooding Attacks[END_REF][START_REF] Wang | SYN-dog: Sniffing SYN Flooding Sources[END_REF], adaptive threshold [START_REF] Bu | Anomaly Network Traffic Detection Based on Auto-Adapted Parameters Method[END_REF], Exponentially Weighted Moving Average (EWMA) [START_REF] Ye | Computer Intrusion Detection Through EWMA for Autocorrelated and Uncorrelated Data[END_REF], Holt-Winters seasonal forecasting [START_REF] Schweller | Reversible Sketches: Enabling Monitoring and Analysis Over High-Speed Data Streams[END_REF], data reduction techniques with Sketch [START_REF] Krishnamurthy | Sketch-based Change Detection: Methods, Evaluation, and Applications[END_REF], SNMP MIB statistical data analysis [START_REF] Yu | Traffic Flooding Attack Detection with SNMP MIB using SVM[END_REF], Principal Component Analysis (PCA) [START_REF] Lakhina | Mining anomalies using traffic feature distributions[END_REF], etc.

Malicious activity usually provokes an abrupt change in the statistical values of the parameters describing the traffic, such as the NetScan produced by worms outbreak, that sends a large number of SYN from the same source IP, to scan the network before propagation phase. In [START_REF] Wang | Change-Point Monitoring for the Detection of DoS Attacks[END_REF], the CUSUM algorithm is used to detect SYN flooding over one time series resulted from aggregating whole traffic in one flow. In [START_REF] Siris | Application of Anomaly Detection Algorithms for Detecting SYN Flooding Attacks[END_REF] a comparison between CUSUM and adaptive threshold for the detection of SYN flooding is presented. Both algorithms have been applied over one time series (whole traffic) for scalability issues and real time processing. Low intensity attacks get smoothed with the normal variations of background traffic in backbone network and pass undetected, i.e. flooding attack with intensity 10 6 packets/s does not produce a noticeable deviation when the total number of normal packets is greater than 10 9 . Furthermore, these methods use static threshold for detecting anomalies, which is not adequate with traffic variations, and may induce false alarm and miss detection. In this paper, we want to overcome these problems (aggregation in one time series and static threshold) through the use of Sketch and dynamic threshold.

The Principal Component Analysis (PCA [START_REF] Lakhina | Mining anomalies using traffic feature distributions[END_REF]) have been proposed and used for dimensionality reduction and anomaly detection. PCA detects anomalies using the evaluation of flows correlations over single link, instead of single time series of whole traffic from many links. However, PCA is very sensitive to the number of dominant PCs. Authors in [START_REF] Khoa | Network Anomaly Detection Using a Commute Distance Based Approach[END_REF] proved that PCA is unable to detect high intensity attacks without the good configuration. In [START_REF] Ringberg | Sensitivity of PCA for Traffic Anomaly Detection[END_REF], the authors showed that methods for tuning PCA are not adequate and starting with a new data set, adjusting parameters is unexpectedly difficult.

Sketch data structure uses the random aggregation for more grained detection than aggregating whole traffic in one time series. It has been used to summarize monitored traffic in a fixed memory, and to provide scalable input for time series analysis. We will exploit the Sketch data structure to derive probability distributions.

Authors in [START_REF] Brauckhoff | Anomaly Extraction in Backbone Networks Using Association Rules[END_REF] experiment the histogram-based detector in order to detect the anomaly behaviors and changes in traffic distributions. They apply Kullback-Leibler divergence between the current and previous measurement distributions. Authors in [START_REF] Tang | Sketch-Based SIP Flooding Detection Using Hellinger Distance[END_REF] apply Hellinger Distance (HD) on Sketch data structure, in order to detect divergence between current and previous distributions of the number of SIP IN-VITE request. In fact, HD must be near zero when probability distributions are similar, and it increases up to one whenever the distributions diverges (e.g. under Invite flooding attacks). In addition, they used the dynamic threshold proposed in [START_REF] Sengar | Detecting VoIP Floods Using the Hellinger Distance[END_REF] during their experimental analysis.

In this paper, we provide a more grained analysis then aggregating the whole traffic in one time series. We develop a general framework that increases the detection accuracy and reduces the false alarm by integrating the Power Divergence over Sketch technique. We show also that KL, HD & χ 2 measures used in [START_REF] Brauckhoff | Anomaly Extraction in Backbone Networks Using Association Rules[END_REF][START_REF] Tang | Sketch-Based SIP Flooding Detection Using Hellinger Distance[END_REF] are special cases of our framework. We will show how PD outperforms these measures by changing its parameter (β).

K-ARY SKETCH

Sketch generates fixed-number of time series [START_REF] Krishnamurthy | Sketch-based Change Detection: Methods, Evaluation, and Applications[END_REF] regardless the number of exiting flows. It allows more grained analysis than aggregating whole traffic in one flow. It is based on random aggregation of traffic in d different hash tables. Let X = x1, x2, . . . , xn denotes the set of input stream, where each item xi = (κi, νi) is identified by a key κi ∈ U , drawn from a fixed universe of items U . νi ∈ R is the value associated with each key. In our model, we use κi = DIP and νi = 1 for each received SYN by the associated Destination IP (DIP).

A Sketch S is a 2D array of d × w cell (as shown in Fig. 1), contains d hash arrays. The arrival of a flow with key κi increments its associated counter in the j th hash table by νi (S j,h j (κ i ) + = νi). The update procedure is realized by d different hash functions, chosen from the set of 2-universal hash functions hj (κi) = {((ajκi + bj) mod PU ) mod w}, to uniformly distribute κi over hash tables and to reduce collisions. The parameter PU is a prime number larger than the maximum number in the universe. aj and bj are random integers smaller than PU , with aj = 0. Using d hash functions, the probability that two keys are aggregated in the same bucket over the d hash tables is (1/w) d . w ( )

, i i v κ i v + i v + i v + d 1 h 2 h ⋮ L h i v + Figure 1: Sketch data structure
Sketch randomly aggregates multiple IP addresses in the same bucket, if the value resulted from hashing the addresses is the same (h(DIP1) = h(DIP2) = j). The counter in each bucket is used to derive a probability as the ratio of counter value to the sum of whole counters in one hash table :

pij = S [i] [j] w j=1 S [i] [j]
(1)

PROPOSED APPROACH

The approach used in this paper to detect flooding attacks is based on Power Divergence. In fact, the idea is to estimate the subjective prior distribution of the traffic and to use it as a baseline probability. Let Pi = [pi,1, ..., pi,w] denotes the probability distribution derived from the i th hash table, in the discrete time interval T . In presence of attacks, the probability distribution changes. One can use this change to detect anomalies. However, with the normal traffic variations, this probability distribution changes also even in the absence of anomaly. This is called false alarms. The objective then is to find a method that detects the attacks and remove the false alarms.

This motivates the need for a quantitative measure of information or more generally a decision theoretic measure of divergence between the basic probability Pi and some other distribution Qi. Power Divergences are generalizations of this decision measure and are associated with strictly convex functions. The Power Divergence has been first defined in [START_REF] Havrda | Quantification method of classification processes: The concept of structural α-entropy[END_REF] and equivalent variants (up to a scale factor β) of this divergence are discussed in [START_REF] Rathie | A Directed-Divergence Function of Type β[END_REF]. The divergence measure is therefore the decision measure that generalizes KL, HD and χ 2 divergence to a broad class of divergence of order β. In fact, the PD is a measure of distance between two probability set of order β given as follows:

P Di = P D(Pi||Qi) = w j=1 pij(pij/qij ) β-1 -1 β(β -1) ( 2
)
where Pi is the probability distribution in the current time interval, and Qi is the probability distribution in previous interval. P Di = 0 iff Pi and Qi are identical (pij = qij ), and P Di > 0 when Pi = Qi. P Di must be near zero under normal traffic, with a large deviation (one spike) when distributions change occurs. d measures (P D1, P D2, . . . , P D d ) are resulted from the d hash tables in Eq. 2. Anomaly induces spike in P Di, and when more than L values of P Di exceed a dynamic threshold, an alarm is raised. The Power Divergence presents some interesting special cases when changing its parameter β, as shown in table 1. Obviously, Power Divergence outperforms then KL, χ 2 and HD measures as they are special cases. In fact, by changing the values of β, one can optimize the detection of attacks compared to the KL, χ 2 and HD measures. By experiments, we will show numerically that for different values of β, the detection efficiency changes. The optimal value of β can then be obtained through our experiments. Pi and Qi in Eq. 1 are derived from the i th hash table in Sketch data structure, in two consecutive discrete intervals. Firstly, the shared counters of the sketch are continuously updated from ongoing traffic during a time interval T . At the end of each interval, the probability pi,j is calculated, and PD is applied to detect deviations.

β Divergence measure -1 1 2 × χ 2 (Qi||Pi) 0 KL(Qi||Pi) 0.5 4 × HD(Pi||Qi) 1 KL(Pi||Qi) 2 1 2 × χ 2 (Pi||Qi)
When L values of P Di are larger than dynamically updated threshold, we raise an alarm. However, PD induces only two spikes (at the start and at the end of attack). As we want to continuously raise alarms for whole duration of the attack, the prior distribution Qi will stop sliding by keeping its value until the end of the attack. However, with the variations of normal traffic, and the similarity of DDoS attacks with flash crowd, we suppose that flooding attacks will span for many intervals to overload a server, in contrast to flash crowd and normal variation. Thus reduce the false alarms. Therefore, we will trigger an alarm only if the deviation lasts more than η intervals.

Detection threshold

Let P D i,kT represent the time series of resulted Power Divergences from different time intervals. To detect deviations in P D i,kT , we derive a subsequent time series P D i,kT that contains the values of P D i,kT smaller than dynamic threshold, i.e. without spikes values. As the empirical rule states that 95% of data fall within 2 standard deviations of mean, P D i,kT contains only values that satisfy:

P D i,kT < µ i,(k-1)T + 2.σ i,(k-1)T (3)
Where µ i,kT & σ i,kT are the mean and the standard deviation of P D i,kT respectively. µ i,kT and σ i,kT are updated dynamically using the EWMA (Exponentially Weighted Moving Average):

µ i,kT = αµ i,(k-1)T + (1 -α)P D i,(k-1)T (4)
σ 2 i,kT = ασ 2 i,(k-1)T + (1 -α)(P D i,kT -µ i,kT ) 2 (5)
The threshold is updated dynamically by adjusting the value of µ i,kT and σ i,kT as shown in eqs. 4 & 5. P D i,kT falls down the threshold (µ i,(k-1)T + 2.σ i,(k-1)T ) under normal condition, and exceeds the dynamic threshold under flooding attack. The decision function for alarms is given in Eq. 6.

To reduce false alarm resulted from normal traffic variations, P D i,kT must exceed the dynamic threshold for η consecutive intervals before raising an alarm.

d(Ai) =      1 if P D i,kT ≥ µ i,(k-1)T + 2.σ i,(k-1)T and η ≥ 3 0 Otherwise (6)

EXPERIMENTAL RESULTS

In this section, we present the performance analysis results for the Power Divergence over Sketch, for SYN flooding attacks detection. Afterward, we conduct analysis to study the impact of parameters on true positive and false alarm ratio. Since web servers use TCP protocol, SYN flooding is the most commonly used attack. Therefore, we focus on presenting the experiment results for SYN flooding detection (due to space limitations), but the proposed approach had applied for the detection of any type of flooding (with TCP, UDP, ICMP, etc.).

We use the real internet MAWI [START_REF]MAWI working group traffic archive[END_REF] trans-Pacific traces from 15/04/2010 12h00 to 18h15 as few hours in the life of The effect of Sketch parameters (width and depth) on the detection accuracy had analyzed in [START_REF] Krishnamurthy | Sketch-based Change Detection: Methods, Evaluation, and Applications[END_REF]. Authors provide a detailed analysis of the impact of the number of hash function (d) and the Sketch width (w) on the detection accuracy. They found that Sketch provides better accuracy when increasing d & w, but at the cost of increasing the required memory and the computational complexity. A trade-off between accuracy and complexity is required, where d = 5 and w = 1024 were chosen as the lowest values to achieve this trade-off in our experiments. Also, the time interval T is subject to reduce the detection delay and the computational complexity. For a value of T = 1sec the computational complexity increases, and for T = 5min the detection delay increases. Therefore, we chose T = 1min as tradeoff.

We analyze these traces using the proposed approach with κi=SIP (Source IP) and vi = 1 for SYN segment. We find many scanning anomalies (SSH Scan, RPC & Netbios Scan, etc.). Afterward, we apply our approach with κi=DIP (Destination IP), and we don't find any existing SYN flooding attacks in these raw traces. Therefore, we inject 9 real DDoS attacks with different intensity to simulate distributed SYN flooding attacks. These attacks are inserted each 30 minutes, and the duration of each one is 10 minutes. The intensity of these attacks begins with a value of 10000 and decreases until 2000 SYN/min. Fig. 2 and Fig. 5 show the variation of the total number of packets before and after the SYN flooding attacks. By comparing these variations, no noticeable deviations between both figures, and inserted attacks don't induce heavy deviations in the time series of the total number of packets. Indeed, Fig. 3 and Fig. 6 show the variation of the number of TCP segments before and after the SYN flooding attacks. We can notice that the shape of traffic variations in both figures is similar. This can be explained by the fact that the intensity of SYN flooding is relatively small when compared to the intensity of the total number of TCP segments. In such cases, the detection of these low intensity attacks is very challenging.

Fig. 4 and Fig. 7 show the variation of SYN before and after SYN flooding attacks. In fact, we can notice the large variations in the total number of SYN (Fig. 4) before attacks injection. Therefore, the aggregation of whole traffic in one time series produces a lot of false alarms for these heavy deviations. On the other hand, the injected SYN flooding attacks produce small variations in the number of SYN (lower than existing normal variations) as shown in Fig. 7. We may not notice the difference between Fig. 4 & Fig. 7 without a deep inspection, given the low intensity of flooding attacks and the large variations in the number of SYN.

We conducted many experiments by varying the values of β in order to find the optimal one. We present the variation 12 with β = 2.5. Whenever the Power Divergence exceeds the dynamic threshold, an alarm is raised. It is important to note the scale difference between these figures, where the intensity of raised pikes increases when increasing the value of β. We don't use the same scale for the sake of clarity.

To evaluate the performance of Power Divergence with different values of β, we investigate the detection rate and the false alarm rate. Receiver Operating Characteristic (ROC) is used for accuracy analysis when varying the value of the threshold. ROC curve shows the variation of the true positive (Eq. 7) in term of false alarm rate (Eq. 8):

DR = T P T P + F N × 100 (7)
Where T P is the number of true positive alerts, and F N is the number of false negative. The false alarm rate is defined as the ratio of false alarms to the number of raised alarms:

F AR = F P T P + F P × 100 [START_REF] Lakhina | Mining anomalies using traffic feature distributions[END_REF] In Fig. 13, we present the ROC curve for different values of β. We found through experiments that for β = 0.5 (PD is proportional to HD), PD achieves a DR=100% with a FAR=30%, for β = 1 (PD is equal to KL), PD achieves a DR=100% with a FAR=23%, for β = 1.5 & β = 2 (PD is proportional to χ 2 ), PD achieves a DR=100% with a FAR=18%, and for β ≥ 2.5, PD achieves a DR=100% with a FAR=10%. It is important to note, that we conduct many experiments with different values of β (with a step 0.5 and up to a value of 100). We found through experiments, with a value of β ≥ 2.5, the intensity of PD for detected attacks increases significantly and proportionally to the attack intensity, but without any change in the ROC. Therefore, PD with β = 2.5 outperforms existing measures and achieves the optimal performance.

CONCLUSION AND PERSPECTIVES

In this paper, we proposed a new framework based on Sketch and Power Divergence for anomaly detection over high speed links. The proposed approach evaluated on real traces with DDoS SYN flooding attacks. Our experimental results showed the capacity of PD in the detection of low intensity attacks. We presented the results of PD with different values of β, as well as the associated ROC curves when varying the threshold. We proved through experiments that PD outperforms existing measures (HD, KL & χ 2 ) when increasing the value of (β). We found that for β = 2.5, PD achieves the optimal performance.

In our future work, we will focus on defense mechanisms, in order to trigger automatic reaction against ongoing attacks. We also intend to provide a method for reducing the amount of monitored data on high speed networks, and to analyze the impact of sampling on the precision of PD divergence measure.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Total number of packets under DDoS attacks

Figure 8 Figure 9 :Figure 10 :

 8910 Figure 8: PD for β = 0.5 (Hellinger Distance)

Figure 11 :

 11 Figure 11: PD for β = 2 (χ 2 divergence)

Figure 12 :

 12 Figure 12: PD for β = 2.5 (Unknown Divergence)

Figure 13 :

 13 Figure 13: Receiver Operating Characteristic (variation of β)

Table 1 :

 1 Special cases of Power Divergence

ACKNOWLEDGMENTS

This research was supported by Korea Science and Engi-neering Foundation, under the World Class University (WCU) program.


Rene Lozi

Ina Taralova

FROM CHAOS TO RANDOMNESS VIA GEOMETRIC UNDERSAMPLING

We propose a new mechanism for undersampling chaotic numbers obtained by the ring coupling of one-dimensional maps. In the case of 2 coupled maps this mechanism allows the building of a PRNG which passes all NIST Test. This new geometric undersampling is very effective for generating 2 parallel streams of pseudo-random numbers, as we show, computing carefully their properties, up to sequences of 12 10 consecutives iterates of the ring coupled mapping which provides more than 10 3.35 10 × × × × random numbers in very short time. Résumé. Nous proposons un nouveau mécanisme de sous-échantillonnage de nombres chaotiques obtenus par couplage en anneau de fonctions unidimensionnelles. Dans le cas de 2 fonctions couplées, ce mécanisme permet de construire un Générateur de Nombres Pseudo-Aléatoires (GNPA) qui satisfait tous les tests NIST. Ce nouveau sous-échantillonnage géométrique est très efficace pour générer deux séries parallèles de nombres pseudo-aléatoires, comme nous le montrons en étudiant très soigneusement leurs propriétés pour des suites de nombres allant jusqu'à 12 10 nombres consécutifs de l'application couplée en anneau, ce qui fournit plus de 10 3.35 10 × × × × nombres aléatoires en un temps très court.

INTRODUCTION

During the last decade, it has been emphasized that the undersampling of sequence of chaotic numbers is an efficient tool in order to build pseudo-random number generators (PRNG) [START_REF] Lozi | Complexity leads to randomness in chaotic systems[END_REF]. Randomness appears to be an emergent property of complex systems of coupled chaotic maps [START_REF] Lozi | Emergence of randomness from chaos[END_REF]. Several kinds of coupling can be considered as ultra-weak coupling or ring coupling,… . An ultra-weak coupling recovers chaotic properties of 1dimensional maps [START_REF] Lozi | Giga-Periodic Orbits for Weakly Coupled Tent and Logistic Discretized Maps[END_REF] when computed with floating numbers or double precision numbers. Chaotic undersampling with thresholds based on one component of the coupled system adds random properties to the chaotic sequences. Double threshold sampled sequence (i.e., using both thresholds of different nature) improves such random properties [START_REF] Lozi | Chaotic Pseudo Random Number Generators via Ultra Weak Coupling of Chaotic Maps and Double Threshold Sampling Sequences[END_REF]. Ring coupling deals when p 1-dimensional maps are constrained on a torus [START_REF] Espinel | New alternate ring-coupled map for multi-random number generation[END_REF], this coupling can directly provide random numbers, without sampling or mixing, provided the number p of maps is large enough, although it is possible to combine these processes with it. However in lower dimension 2 and 3, the chaotic numbers are not equidistributed on the torus. Therefore we introduce a particular "geometric" undersampling based on the property of piecewise linearity of the invariant measure of system. This new geometric undersampling is very effective for generating parallel streams of pseudo-random numbers with a very compact mapping.

In Section 1 we briefly recall properties of chaotic mappings, when used alone or ultraweakly coupled. Section 2 describes the route from chaos to randomness via chaotic undersampling, discovered during the last decade. In Section 3, we introduce geometric undersampling in the scope of ring coupled mapping.

COMPUTATION OF CHAOTIC NUMBERS

Disappointing chaotic numbers

Chaos theory studies the behavior of dynamical systems that are highly sensitive to initial conditions, an effect which is popularly referred to as the butterfly effect. Small differences in initial conditions (such as those due to rounding errors in numerical computation) yield widely diverging outcomes for chaotic systems, rendering long-term prediction impossible in general. This happens even though these systems are deterministic, meaning that their future behavior is fully determined by their initial conditions, with no random elements involved. In other words, the deterministic nature of these systems does not make them predictable. The first example of such chaotic continuous system in the dissipative case was pointed out by the meteorologist E. Lorenz in 1963 [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF].

In order to study numerically the properties of the Lorenz attractor, M. Hénon an astronomer of the observatory of Nice, France, introduced in 1976 a simplified model of the Poincaré map of this attractor [START_REF] Hénon | A Two-dimensional mapping with a strange attractor[END_REF]. The Lorenz attractor being imbedded in dimension 3, the corresponding Poincaré map is a mapping from the plane 1² into 1². Hence the Hénon mapping is also defined in dimension 2 and is associated to the dynamical system x y ax y bx

+ + 1 = + - 2 = 3 , (1)
which has been extensively studied since thirty six years. More simple dynamical systems in dimension one, on the interval

[ ] 1,1 J = - ⊂ 1 into itself 1 n a n x f ( x ) + = , (2)
f L ( x ) ax ≡ = - , (3)
or the symmetric tent map

1 a a f T ( x ) a x ≡ = - , (4)
have also been fully explored in the hope of generating random numbers easily [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF]. However when a dynamical system is realized on a computer using floating point or double precision numbers, the computation is of a discretization, where finite machine arithmetic replaces continuum state space. For chaotic dynamical systems in small dimension, the discretization often has collapsing effects to a fixed point or to short cycles [START_REF] Lanford | Some informal remarks on the orbit structure of discrete approximations to chaotic maps[END_REF][START_REF] Gora | Absolutely continuous invariant measures that cannot be observed experimentally[END_REF].

It seems that the computation of numerical approximations of the periodic orbits leads to unpredictable and somewhat enigmatic results. As says O. E. Lanford III [START_REF] Lanford | Some informal remarks on the orbit structure of discrete approximations to chaotic maps[END_REF] "The reason is that because of the expansivity of the mapping the growth of roundoff error normally means that the computed orbit will remain near the true orbit with the chosen initial condition only for a relatively small number of steps typically of the order of the number of bits of precision with which the calculation is done. It is true that the above mapping like many 'chaotic' mappings satisfies a shadowing theorem (see [START_REF] Palmer | Shadowing in Dynamical Systems: Theory and Applications[END_REF][START_REF] Pilyugin | Shadowing in Dynamical Systems[END_REF]) which ensures that the computed orbit stays near to some true orbit over arbitrarily large numbers of steps. The flaw in this idea as an explanation of the behavior of computed orbits is that the shadowing theorem says that the computed orbit approximates some true orbit but not necessarily that it approximates a typical one."

The collapsing of iterates of dynamical systems or at least the existence of very short periodic orbits, their non constant invariant measure, and the easily recognized shape of the function in the phase space avoid the use of one-dimensional map (logistic, baker, or tent, …) as a Pseudo Random Number Generator (see [START_REF] Lozi | Can we trust in numerical computations of chaotic solutions of dynamical systems[END_REF] for a survey).

Remark 1.1 However, the very simple implementation in computer program of chaotic dynamical systems led some authors to use it as a base of cryptosystem [START_REF] Baptista | Cryptography with chaos[END_REF][START_REF] Ariffin | Modified Baptista type chaotic cryptosystem via matrix secret key[END_REF]. In addition it seems that for some applications, chaotic numbers are more efficient than random numbers, that is the case for evolutionary algorithms [START_REF] Pluhacek | Extended Initial Study on the Performance of Enhanced PSO Algorithm with Lozi Chaotic Map[END_REF][START_REF] Tang | Parrondo's games with chaotic switching[END_REF] or chaotic optimization [START_REF] Araujo | Particle swarm approaches using Lozi map chaotic sequences to fuzzy modelling of an experimental thermal-vacuum system[END_REF].

In this paper we show how to overcome the poor quality of chaotic generators using geometric undersampling. This special undersampling we introduce in this article is one among others undersampling processes we have studied before. In order to explain the difference between these processes we give in Sec. 2 a brief survey of them. Before doing this survey, we have to show how to stabilize the chaotic properties of chaotic number when realized on a computer. 1.2.

Long periodic orbits for ultra-weakly coupled tent map

The first step in order to preserve the genuine chaotic properties of the continuous models in numerical experiments is reached considering ultra-weak multidimensional coupling of p one-dimensional dynamical systems instead of solely a one-dimensional map.

System of 2-coupled symmetric tent map

In order to simplify the presentation below, we use as an example the symmetric tent map (4) with the parameter value a = 2, later denoted simply as f, even though others chaotic map of the interval (as the logistic map, the baker transform, …) can be used for the same purpose (as a matter of course, the invariant measure of the chaotic map considered is preserved).

When p = 2, the system is simply described by Eq. ( 5)

1 1 1 1 2 2 1 1 n n n n n n x ( 1 ) f ( x ) 1 f ( y ) y 1 f ( x ) ( 1 ) f ( y ) + + = - + 1 2 = + - 3 , (5)
We use generally for double precision numbers. In both cases, with these numerical values, the collapsing effect disappears and the invariant measure of any component is the Lebesgue measure [START_REF] Lozi | Giga-Periodic Orbits for Weakly Coupled Tent and Logistic Discretized Maps[END_REF] as we show below. In the case of computation using floating points, starting form most initial condition, it is possible to find a Mega-Periodic orbit (i.e. with period equal to 1,320,752). When computations are done with double precision number it is not possible to find any periodic orbit, up to 11 [START_REF] Gora | Absolutely continuous invariant measures that cannot be observed experimentally[END_REF] 10 n = × iterations. In [START_REF] Lozi | Giga-Periodic Orbits for Weakly Coupled Tent and Logistic Discretized Maps[END_REF] the computations have been performed on a Dell computer with a Pentium IV microprocessor using a Borland C compiler computing with ordinary (IEEE-754) double precision numbers. ε converges towards 0, the iterates of each component x n and y n of equation [START_REF] Gora | Absolutely continuous invariant measures that cannot be observed experimentally[END_REF] converge to the Lebesgue measure (Fig. 1).

System of p-coupled symmetric tent map

More generally, the coupling of p maps takes the form ( ) ( )

1 n n n X F X A f ( X ) + = = ⋅ , (6) where 1 ( ) ( ) (
=1-1 1 1 1 1 =1- 1 1 1 A 1 1 =1- 1 ε ε ε 3 3 2 4 2 2 2 4 2 2
3 3 = = = = ≠ = - = 4
ε ε = = ≠
A on the diagonal (the matrix A is always a stochastic matrix iff the coupling constants verify 0 i , j ε > for every i and j).

It is noteworthy that these families of very weakly coupled maps are more powerful than the usual formulas used to generate chaotic sequences, mainly because only additions and multiplications are used in the computation process, no division being required. Moreover the computations are done using floating point or double precision numbers, allowing the use of the powerful Floating Point Unit (FPU) of the modern microprocessors. In addition, a large part of the computations can be parallelized taking advantage of the multicore microprocessors which appear on the market of laptop computers.

Moreover, a determining property of such coupled map is the high number of parameters used (

( 1) p p ×for p coupled equations) which allows to choose it as cipher-keys, when used in chaos based cryptographic algorithms, due to the high sensitivity to the parameters values [START_REF] Lozi | Emergence of randomness from chaos[END_REF]. It at be shown also, that using control theory techniques, synchronization of two systems [START_REF] Hénaff | Dynamical analysis of a new statistically highly performant deterministic function for chaotic signals generation[END_REF], with p = 2 or 3, can be reached via exact (dead-beat) or asymptotic observers [START_REF] Hénaff | Exact and Asymptotic Synchronization of a new weakly coupled maps system[END_REF]. 1.2.3.

Approximated distribution function

In order to assess numerical computations more accurately, we define an approximation , ( )

M N P
x of the invariant measure also called the probability distribution function linked to the 1-dimensional map f, when computed with floating numbers (or numbers in double precision). In this scope we consider a regular partition of M small intervals (boxes) r i of J defined by Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

6 i 2i s 1 , i 0, M M = -+ = (9) i i i 1 r s , s , i 0, M -2 + = = B B C C and [ ] M 1 M 1 r s , 1 - - = (10)
the length of each box is equal to 2 M and the i r intervals form a partition of the interval (11) All iterates f (n) (x) belonging to these boxes are collected, after a transient regime of Q iterations decided a priori, (i.e. the first Q iterates are neglected). Once the computation of N + Q iterates is completed, the relative number of iterates with respect to N/M in each box r i represents the value ( )

J 1 0 M i J r - = 5
N i
P s . The approximated ( ) N P x defined in this article is then a step function, with M steps. As M may vary, we define ( )

, ( ) # M N i i M P s r N = ( 12
)
where #r i is the number of iterates belonging to the interval r i . ) (

, x P N M is normalized to 2 on the interval J . i i N M N M r x s P x P ∈ ∀ = ) ( ) ( , , (13)
In the case of p-coupled maps, we are more interested by the distribution of each component ( ) , , , ,

1 2 1 p 2 x x x x 6
of X rather than the distribution of the variable X itself in p J .

We then consider the approximated probability distribution function , ( )

j M N P
x associated to one among several components of F(X) defined by [START_REF] Hénaff | Dynamical analysis of a new statistically highly performant deterministic function for chaotic signals generation[END_REF], which are one-dimensional maps. In this paper we use equally disc N for M and iter N for N, when they are more explicit.

The discrepancies 1

E (in norm 1 L ), 2 E (in norm 2 L
) and E ∞ (in norm L ∞ ) between disc iter j N , N P ( x ) and the Lebesgue measure, which is the invariant measure associated to the symmetric tent map, are defined by

1 1, , , ( ) ( ) 1
disc iter disc iter j j N N N N L E x P x = - (14) 2 2, , , ( ) ( )
1 disc iter disc iter j j N N N N L E x P x = - (15) , , , ( ) ( )
1 disc iter disc iter j j N N N N L E x P x ∞ ∞ = - (16)
As previously said, Fig. 1 shows the convergence of the density of iterates of the components of 2-coupled symmetric tent maps to the Lebesgue measure when 1 ε converges towards 0. Moreover, for a fixed value of disc N when the number iter N increases, the discrepancy between disc iter j N , N P ( x ) and the Lebesgue measure is expected to converge towards 0, except if there exist periodic orbits of finite length lower than iter N which capture the iterates. In this case whatsoever the value of iter N is, the approximated distribution function converges to the distribution function of the periodic orbit, if it is unique, or to the average of the distribution functions of the periodic orbits observed, if not. ( )

disc iter N N E
x versus the number of iterates of the approximated distribution functions, with respect to the first variable 1 x , for 2 and 3-coupled symmetric tent map. Same results are obtained for the other variables 2 x or 3

x . The 3-coupled symmetric tent maps model considered here with very very small value of ε 1 , seems a sterling model of generator of chaotic numbers with a uniform distribution of these numbers over the interval J . It produces very long periodic orbits: Gigaperiodic orbits (i.e. with length of period between 9 10 and 12 10 ) when computed with simple precision numbers, and orbits of unknown length when computed with double precision numbers. However these chaotic sequences are not at all random sequences.

THE ROUTE FROM CHAOS TO RANDOMNESS VIA CHAOTIC UNDERSAMPLING

Chaotic numbers are not pseudo-random numbers because the plot of the couples of any component ( )

l n l n x x 1 , + of iterated points ( ) 1 ,
n n X X + in the corresponding phase plane reveals the map f used as one-dimensional dynamical systems to generate them via Eq. ( 6). Nevertheless, we have recently introduced a family of enhanced Chaotic Pseudo Random Number Generators (CPRNG) in order to compute faster long series of pseudorandom numbers with desktop computer [START_REF] Lozi | Chaotic Pseudo Random Number Generators via Ultra Weak Coupling of Chaotic Maps and Double Threshold Sampling Sequences[END_REF][START_REF] Lozi | Complexity leads to randomness in chaotic systems[END_REF]. This family is based on the previous ultra weak coupling which is improved in order to conceal the chaotic genuine function.

In this Section, we describe briefly how works this first process of undersampling, the chaotic one.

Chaotic undersampling

In order to hide f in the phase space ( )

l l n n 1
x , x + two mechanisms are used. The pivotal idea of the first one mechanism is to sample chaotically the sequence ( ) x , is strictly greater (or smaller) than a threshold T J ∈ , with l ≠ m, for 1 ≤ l, m ≤ p.

That is to say to extract the subsequence ( )

( 0 ) (1) ( 2) ( ) ( 1)
, , , , , ,

q q l l l l l n n n n n x x x x x + 6
6 denoted here ( )

3 3 , , , , , , 1 2 1 0 + q q x x x x x
of the original one, in the following way

given 1 , , l m p l m ≤ ≤ ≠ { } ( ) ( 1)
( )

,

q l m q n q q r r n x x with n Min r n x T - - ∈ = - 1 D 2 = = > > D 3 7 (17)
The sequence ( )

3 3 , , , , , , 1 2 1 0 + q q x x x x x
is then the sequence of chaotic pseudo-random numbers.

The above mathematical formula can be best understood in algorithmic way. The pseudocode, for computing iterates of ( 17) corresponding to N iterates of ( 6) is:

( )

1 2 p 1 p 0 0 0 0 0 X x , x , , x , x seed - = = 6 n 0; q 0 = = do { while n < N do { while ( ) m n
x T ≤ compute ( )

1 2 p 1 p n n n n
x , x , ,

; n + +} compute ( )

1 2 p 1 p n n n n
x , x , , x , x -6 ; then

( ) q l q n
n( q ) n ; x x ; = = n + +; q + +} This chaotic sampling is possible due to the independence of each component of the iterated points n X versus the others [START_REF] Lozi | New Enhanced Chaotic Number Generators[END_REF].

Remark 2.1. Albeit the number iter NSampl of pseudo-random numbers q x corresponding to the computation of N iterates is not known a priori, considering that the selecting process is again linked to the uniform distribution of the iterates of the tent map on J , this number is equivalent to 2 1

N T - .

Chaotic mixing

A second mechanism can improve the unpredictability of the pseudo-random sequence generated as above, using synergistically all the components of the vector n X , instead of two. Given p -1 thresholds

1 2 1 p T T T J - < < < ∈ 3 ( 18
)
and the corresponding partition (note that this partition of J is different from the regular previous one [START_REF] Lozi | Giga-Periodic Orbits for Weakly Coupled Tent and Logistic Discretized Maps[END_REF] used for the approximated distribution function).

[ ] 0 1 J 1,T = - , ] [ 1 1 2 J T ,T = , [ [ 1 1 1 k k k J T ,T for k p + = < < -and 1 1 1 p p J T , - - - - - - - - B B B B B B B B = = = = C C C C C C C C , with 1 0 p k k J J - = = 5 (19
The second simple mechanism is based on the chaotic undersampling combined with a chaotic mixing of the p-1 sequences ( )

6 6
, , , , , ,

1 1 1 1 2 1 1 1 0 + n n x x x x x , ( ) 6 6 , , , , , , 2 1 2 2 2 2 1 2 0 + n n x x x x x ,…, ( ) 6 6
, , , , , ,

1 1 1 1 2 1 1 1 0 - + - - - - p n p n p p p x x x x x ,….
using the last one ( )

6 6 , , , , , , 1 2 1 0 p n p n p p p x x x x x
+ in order to distribute the iterated points with respect to this given partition, defining the subsequence ( )

3 3 , , , , , , 1 2 1 0 + q q x x x x x
(in pseudocode) by ( )

1 2 p 1 p 0 0 0 0 0 X x , x , , x , x seed - = = 6 n 0; q 0 = = do { while n < N do {while ( ) p n 0 x J ∈ compute ( ) 1 2 p 1 p n n n n
x , x , , x , x -

6

; n + +} compute ( )

1 2 p 1 p n n n n x , x , , x , x - 6 let k be such that p n k x J ∈ ; then ( ) q k q n
n( q ) n ; x x = = ; n + +; q + +} Remark 2.2. In this case also, iter NSampl is not known a priori, however, considering that the selecting process is linked to the uniform distribution of the iterates of the tent map on J , one has Remark 2.3. This second mechanism is more or less linked to the whitening process [START_REF] Viega | Practical random number generation in software[END_REF][START_REF] Viega | Secure programming cook book for C and C++[END_REF].

Remark 2.4. Actually, one can choose any of the components in order to sample and mix the sequence, not only the last one.

Enhanced chaotic undersampling

On can eventually improve the CPRG previously introduced with respect to the infinity norm instead of the 1 L or 2 L norms because the L ∞ norm is more sensitive than the others ones to reveal the concealed f [START_REF] Lozi | Chaotic Pseudo Random Number Generators via Ultra Weak Coupling of Chaotic Maps and Double Threshold Sampling Sequences[END_REF]. For this purpose we introduce a second kind of threshold ' T ∈ 7 , together with 1 1 , , p T T J -∈ 3 such that the subsequence ( ) n( q ) n ; x x = = ; n + +; q + +} Remark 2.5. In this case also, iter NSampl is not known a priori, it is more complicated to give an equivalent to it. However, considering that the selecting process is linked to the uniform distribution of the iterates of the tent map on J , and to the second threshold T', it comes that 1 2 ,

3 3 , , , , , , 1 2 1 0 + q q x x x x x is defined (in pseudo-code) by ( ) 1 2 p 1 p 0 0 0 0 0 X x , x , , x , x seed - = = 6 n 0, q 0 = = do { while n < N do {while ( ) p ( q 1 ) n 0 n n T ' and x J - ≤ + ∈ compute ( ) 1
1 iter N N NSampl Min T T 1 E ≤ 2 F ′ - 3 .
Remark 2.6. The second kind of threshold ' T can also be used with only the chaotic sampling, without the chaotic mixing.

A window of emergence of randomness

In [START_REF] Lozi | Complexity leads to randomness in chaotic systems[END_REF][START_REF] Lozi | Emergence of randomness from chaos[END_REF] we show that if one consider the errors

1 1, , , ( ) ( )
1 disc iter disc iter N N N N L E x P x = - , 2 2, , , ( ) ( ) 1
disc iter disc iter N N N N L E x P x = - and , , , ( ) ( ) 1
- - - - - - - - B B B B ∈ ∈ ∈ ∈ C C C C
, in the case p = 4 and

1 i , j i i ε ε ε = = .
We have also performed NIST test developed by the National Institute of Standards and Technology [START_REF] Rukhin | A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications[END_REF], in order to check carefully the random nature of such numbers [START_REF] Hénaff | Dynamical analysis of a new statistically highly performant deterministic function for chaotic signals generation[END_REF].

Then there is a route from chaos to randomness using the process of chaotic undersampling.

GEOMETRIC UNDERSAMPLING

The previous route from chaos to randomness uses chaotic undersampling. It is possible to couple in another way p tent maps on the torus [ ] 1,1 p p p J = -⊂ 1 , which can directly provide random numbers, without sampling or mixing, provided p is large enough, although it is possible to combine these processes with it. After reviewing this ring coupling in high dimension, we introduce the new geometric undersampling in order to obtain randomness with small values of p (for example p = 2).

Ring coupled mapping

Consider the mapping defined on the p-dimensional torus :

p p p M J J → 2 1 1 1 n 1 n n n 1 3 2 2 2 n 2 n n 1 n p p p 1 m n n 1 n m n k x 1 2 x x x k x 1 2 x x x M x x k x 1
+ + + + + + + + 1 1 1 1 < - < - < - < - 2 2 2 2 > > > > 3 3 3 3 (21)
The particularity of this coupling is that each variable j x is coupled only with itself and j 1

x + + + + , as displayed on Fig. 3a. At first glance, in order to enrich the random properties of the map, it could seem interesting to add supplementary cross couplings between these variables, as shown on Fig. 3b. However in this case a cross-coupling is inappropriate because it would increase the determinism and therefore deteriorate the statistical properties which we are looking for.

To evaluate the random properties of these generators, the set of NIST tests have been used again.

Figure 3a (left). Ring coupling between the variables j

x .

Figure 3b (right). Cross coupling between the variables j

x .

The random properties validations of both a 4-dimensional system and a 10-dimensional one have been carried out [START_REF] Espinel | New alternate ring-coupled map for multi-random number generation[END_REF]. For this purpose, the chaotic carrier output needs to be quantised and binarized (0 and 1) in order to be validated as being random using NIST tests. Therefore, different methods of binarization (converting real signals to binary ones) have been implemented and compared.

A first 1-bit binarization has been applied to the system (21) output, defined as

j n n y x = = = = with 8 9 1, j p ∈ n if ( y 0 ) b 1 else b 0 ≥ = ≥ = ≥ = ≥ = 1 1 1 1 2 2 2 2 = = = = 3 3 3 3 ( 22
)
The results showed to be highly sensitive to the type of binarization. Eventually, after testing several different methods, a 32-bit binarization has been chosen as being the most suitable solution. Because the system is confined to the p-dimensional torus p J , 31 bits are assigned to represent the decimal part, and 1 bit to the sign. To illustrate the results, the NIST tests for the 4-dimensional system with parameters 10 bits, length of bit string = 6 10 bits, quantity of bit strings = 100. The output of the system has been arbitrary chosen as being:

4 n y x = .
Furthermore, as the results show their independence from the initial conditions, every bit string in this test is the resulting sequence of a different randomly chosen initial condition. The criterion for a successful test is that the p-value has to be superior to the significance level (0.01 for this case). For the present model, all tests were successful thus the sequences can be accepted as being random.

+ + + + = - = = - = = - = = - =
, each sequence of components satisfies the NIST test for randomness.

Low dimensional model

Although the system ( 20) is a good PRNG when 4 p ≥ , in lower dimension 2 and 3, the chaotic numbers are not equidistributed on the torus (see Fig. 5).

In order to improve the ring coupling mechanism in low dimension, we introduce now a very new type of undersampling based on geometric nature of the invariant measure. We present this very new mechanism which allows the emergence of randomness from chaos, in the simplest case, the 2-dimensional ring mapping 2 M on the square 2 J , with

1 2 k k 1 = = = = = = = = .
Let 2 M be defined by

1 1 2 n 1 n n 2 2 1 n 1 n n x 1 2 x x x 1 2 x x + + + + + + + + 1 1 1 1 = - + = - + = - + = - + D D D D 2 2 2 2 = - + = - + = - + = - + D D D D 3 3 3 3 (23) with j n 1 j n 1 if ( x 1 ) add 2 if ( x 1 ) substract 2 + + + + + + + + 1 1 1 1 < - < - < - < - 2 2 2 2 > > > > 3 3 3 3 (24)
3.2.1. Critical lines.

Figure 5 shows the distribution of the iterates of system (23) (the transient of the first 6 10 iterations has been cut off). It can be observed that the attractor contains regions where the point density is lower, and two lozenge-like holes. It is possible to define critical lines which form a partition of the square 2 J . The critical lines CL [START_REF] Mira | Chaotic dynamics in two-dimensional noninvertible maps[END_REF] are singularities of dimension 1 and represent an important tool for the analysis of noninvertible maps. The holes on Fig. 5 are completely delimited by segments of the critical lines The critical lines separate regions of the phase space with different number of preimages (backward iterates). In the case of piecewise linear maps, they are the first iterates of the lines of discontinuity 1 CLof the system. For the two dimensional system [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF], there are four groups of critical lines CL with preimages

Critical lines A:

1 1 0 A CL : x - = and 2 1 2 1 1 2 2 1 2 1 2 1 0 2 1 0 A A if : x x x CL : if CL x x x 1 = - - > D 2 = - < D 3 (25)
Critical lines B:

1 1 1 B CL : x - = - and [ ] [ ] [ ] [ ] 1 2 1 2 1 1 2 2 1 2 1 1 3 2 1 2 1 1 4 2 1 2 1 1 2 0 0 0 5 2 2 0 1 0 5 2 2 0 0 5 1 2 0 0 5 0 B B B B : if , , . CL x x x x : if , , . CL x x x x : if , . , CL x x x x : if , . , CL x x x x 1 = < ∈ D = - - > ∈ -- D 2 = - < ∈ D D = - > ∈ - 3 (26)
Critical lines C: J (a square) [START_REF] Taralova | Dynamical and statistical analysis of a new Lozi function for random numbers generation[END_REF].

2 1 0 C CL : x - = and ( ) ( ) 1 2 1 1 1 2 2 1 1 1 1 1 0 2 1 1 0 2 C C : if CL x x x : if CL x x x 1 = - + > D D 2 D = + < D 3 (27)

Critical lines D:

Markov partition.

Our aim is first to use the partition defined by these critical lines in order to do a cell-tocell analysis and, by the means of a Markov process, to compute explicitly the invariant measure of iterates associated to system [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF]. Figure 6 displays the 32 sub-regions of the square 2 J , labelled from a to p and a' to p'. For clarity of the presentation, we have labelled from (I) to (IV), the four quadrants of 2 J . 15

Straightforward computation shows that the images of each region, by the mapping 2 M , is either one, two or three regions of the same partition of the square 2 J . Figures 7a and 7b display the images of the regions embedded in the first quadrant (I). Figures 8a and8b display the images of the regions embedded in the second quadrant (II). The colour is the same for every region and its corresponding image, except when two regions are mapped on the same region, in this case there is a mix of colours on the common part of the image. The overall correspondence between regions of the partition and their image is given by the Markov matrix a M which is displayed on Table 1. The computation of the coefficients of this matrix, which are rational numbers, is based on the ratios of surfaces of bounded regions.

1 = = = = D D D = = = = 2 D D = = = = D 3 (29)

Exact computation of invariant measure associated to 2 M

With the help of Markov matrix a M , it is straightforward to compute explicitly the invariant measure associated to 2 M . For every region on Figure 6, we define a quantity of initial points called i Q , i 1,32 = = = = uniformly scattered on it, and we compute its surface i S .

We normalise both quantities to

i i Q Q 4 = = = = = = = = A A A A
, and

i i S S 4 = = = = = = = = A A A A
. Hence it is possible to define the density of iterates on each region. the vectors of quantities and densities obtained applying (30) to every region. Then starting from an arbitrary initial repartition of points on 2 J , say

i i i Q d S = = = = (30) Let 1 32 Q Q Q 2 4 5 4
1 0 0 32 0 Q Q Q 4 5 4
Q M Q + + + + = = = = (31)
The sequence of vectors {

{ { { } } } } m m Q 7 ∈ ∈ ∈ ∈ converges to a limit vector Q which satisfies t a Q M Q = = = = (32)
and gives the invariant measure, the density of which is the vector D , using (30).

Numerical results:

Starting from

1 0 0 32 0 Q 1 / 8 Q Q 1 / 8

Geometric undersampling

The exact computation of the density D of the invariant measure shows that this density is constant on each region. The geometric undersampling process consists in magnifying a square G included in one region (as for example the square

Algorithm of geometric undersampling

Let

1 1 2 2 l r l r G x , x x , x B B = × C C
the square in which we will undersample the iterate of ( 23) and, NSampl the number of geometrically undersampled iterates is not known a priori, however, considering that the selecting process is linked to the uniform distribution of the iterates of the tent map on 2 J , one has ( )

Numerical tests

We have applied this process in the case of the square G of Fig. 9, with J . We have also used NIST test to confirm the random property of the geometrical undersampling process. They are all successful (Fig. 12).

CONCLUSION

We have proposed a new mechanism of undersampling of chaotic number obtained by the the ring coupling mechanism of one-dimensional maps. In the case of 2 coupled maps this mechanism allows the building of a PRNG which passes all NIST Test.

This new geometric undersampling is very effective for generating 2 parallel streams of pseudo-random numbers, as we have shown, computing carefully their properties, up to sequences of 12 10 consecutives iterates of ( 23) which provides more than

Figure 1 .

 1 Figure 1. Density of iterates of 2-coupled symmetric tent maps, double precision, 5 10 disc N = , 2

Figure 2

 2 Figure 2 shows the errors 1 1, ,

  Figure 2. Error

+

  generated by the l-th component l x , selecting l n x every time the value m n x of the m-th component m

  correlated distribution functions which assess the independence of each component of the iterated points n X versus the others, a window of emergence comes clearly into sight for the values

1 (

 1 The chosen conditions are: length of the original sequence =8

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

Figure 4 .

 4 Figure 4. Example of NIST Test for i

1 CL

 1 given by Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

14 Figure 5 . 2 M

 1452 Figure 5. Critical lines of the map

Figure 6 .

 6 Figure 6. The 32 sub-regions for a partition of the square 2J .

Figure 7 a

 7 Figure 7 a (left). The nine regions a to i of quadrant (I) Figure 7 b (right). The images 2 M (I) of the nine regions of quadrant (I).

Figure 8 a

 8 Figure 8 a (left). The seven regions j to p of quadrant (II) Figure 8 b (right). The images 2 M (II) of the seven regions of quadrant (II).

  36,0.64 0.36,0.64 = × included in region m on Fig.9), up to the size of the square 2 J .

Figure 9

 9 Figure 9. The square

  algorithmic form, the pseudo-code, to geometric n < N compute ( )

d

  is the density of the measure in region m.

Figure 10 a

 10 Figure 10 a (left). densities of the seven regions j, k, l, m, n, o, p of quadrant (II) Figure 10 b (right). Uniform density of iterates in the square [ ] [ ] G 0.36,0.64 0.36,0.64 = × of quadrant (II).

Figure 12 .

 12 Figure 12. Geometrical undersampling: each sequence of components satisfies the NIST test for randomness.

  short time. In a forthcoming paper we will test the 3-dimensional case.

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

	corresponding to the logistic map
	a	a	1	2
				3

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

  )

			f X	n	4 6 6 = 6 6 6 8	n p n f x f x 2 2	5 7 7 7 7 7 9	,	X	n	1 n p n x x 4 5 6 7 6 7 = 6 7 8 9 6 7 6 7 2 2	, (7)
	and
	j p
	1,1	1,	j		1,2					1,p-1	1,p
	j	2
						j p
	2,1		2,2			2,	j			2,p-1	2,p
					j	1, 2 j
											j p	1
	p,1										p,p-1	, p p	, p j
											j	1

)

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

	then	q		k ( ) q n
	compute (	2 n x , x , , x , x p 1 n n -6 n p ) 2 p 1 p n n n x , x , , x , x 1 n -6	; n + +}
	let k be such that p n x	∈	J	k
					10

Table 1 .

 1 Markov Matrix a M .

	In order to display the 32 32 × × × ×	matrix a M on one page, we have labelled the coefficients
	using letters which are not related to the names of the regions.
	o	1 12	; p		1 20	;q	3 20	;r	4 20	;
		s	1 9	;t		2 9	;u	4 9	;v	1 6	;
		w	1 3	; x	1 5	; y	3 5	; z	2 3	;

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).

	500 Q	= = = =	4 4 4 4 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 8 8 8 8 6 6 6 6 6 6	1 500 2 500 3 500 4 500 29 500 30 500 31 500 32 500 Q Q Q Q Q Q Q Q 2	= = = =	4 4 4 4 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 8 8 8 8	1 / 14 3 / 28 1 / 14 3 / 28 4 / 7 3 / 28 3 / 28 1 / 7 2	5 5 5 5 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 9 9 9 9	= = = =	Q	, which gives using (30),	500 D	d d d d	1 500 2 500 3 500 500 32 31 500 500 30 500 29 4 500	10 / 7 5 / 7 10 / 7 5 / 7	D .
	Remark 3.1. Computing directly this density, iterating (23) up to 10 11 iterates leads to the
	same result.
							= = = =	4 4 4 4 6 6 6 6 6 6 6 6 6 6 6 6 8 8 8 8	2	5 5 5 5 4 4 4 4 7 7 7 7 6 6 6 6 = = = = 7 7 7 7 6 6 6 6 6 6 6 6 7 7 7 7 8 8 8 8 9 9 9 9	2	5 5 5 5 7 7 7 7 7 7 7 7 7 7 7 9 9 9 9 7	, Q is obtained rapidly, as

  In this case, the undersampling process provides two streams of pseudo-random numbers. Remark 3.3. In this case,

					1 n x , x ;if ( ) 2 n 1 2 n n x , x	G ∈ then
	1 q x 2 =	B C	1 n x x mean 1 1 1 r l x x --	2 q , x 2 =	B C	2 n x x mean 2 2 2 r l x x --	; q=q+1; n + +}
	Remark 3.2. iter

  Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings). II). On Fig.11the square is magnified up to the size of the square2 J . The vertical scale is fitted near the invariant Lebesgue measure.

	Figure 10b shows the uniform density of iterates in the square [ ] [ ] G 0.36,0.64 0.36,0.64 = × 25 29 33 37 41 45 49 3200 3250 3300 3350 3400 3450 3500 3550 of quadrant (1 3600
	5
	9														21
	13	17	21	25	29	33	37	41	45	49	1	5	9	13	17
	3200-3250	3250-3300	3300-3350	3350-3400	3400-3450	3450-3500		3500-3550	3550-3600
	Figure 11. Uniform density of iterates of the square [ ] [ ] G 0.36,0.64 0.36,0.64 = × magnified to the square 2

Université de Nice Sophia-Antipolis, Laboratoire J. A. Dieudonné, UMR CNRS 7351, Parc Valrose, 06108 NICE, Cedex 02, France, email : rlozi@unice.fr

L'UNAM, IRCCyN, UMR CNRS 6597, Ecole Centrale de Nantes, 1, rue de la Noë, BP 92101, 44321, NANTES Cedex

3, France, email : Ina.Taralova@irccyn.ec-nantes.fr

Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings). Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings). Submitted for publication In Proceedings of European Conference on Iteration Theory (ECIT 2012), Ponta Delgada, Açores, Portugal, September 9-15, 2012. (to be published in European Series in Applied and Indsutrial Mathematics, ESAIM : Proceedings).


Qiang Bai

Romain Desprat

Bernard Klein

Jean-Marc Lemaitre

John De Vos
email: <john.devos@inserm.fr

Embryonic stem cells or induced pluripotent stem cells? A DNA integrity perspective

Keywords: Animals, Cell Culture Techniques, Cell Differentiation, DNA, genetics, Embryonic Stem Cells, cytology, immunology, Epigenesis, Genetic, Genomic Instability, Humans, Induced Pluripotent Stem Cells, cytology, immunology, Mice Induced pluripotent stem cells, cell reprogramming, genetic abnormalities, pluripotency, DNA damage, genome integrity

Induced pluripotent stem cells (iPSCs) and embryonic stem cells (ESCs

) are two types of pluripotent stem cells that hold great promise for biomedical research and medical applications. iPSCs were initially favorably compared to ESCs. This view was first based on ethical arguments (the generation of iPSCs does not require the destruction of an embryo) and on immunological reasons (it is easier to derive patient HLA-matched iPSCs than ESCs). However, several reports suggest that iPSCs might be characterized by higher occurrence of epigenetic and genetic aberrations than ESCs as a consequence of the reprogramming process. We focus here on the DNA integrity of pluripotent stem cells and examine the three main sources of genomic abnormalities in iPSCs: (1) genomic variety of the parental cells, (2) cell reprogramming, and (3) cell culture. Recent reports claim that it is possible to generate in vitro mouse or human iPSC lines with a mutation level similar to that of the parental cells, suggesting that genome-friendly " " reprogramming techniques can be developed. The issue of iPSC DNA integrity clearly highlights the crucial need of guidelines to define the acceptable level of genomic integrity of pluripotent stem cells for biomedical applications. We discuss here the main issues that such guidelines should address.

INTRODUCTION

Embryonic stem cells (ESCs), which can theoretically differentiate into any cell type, have become a key cellular tool in biology and have paved the way for cell-based treatments of human diseases 1 . However, human ES cells (hESCs) have two major drawbacks. The [ ] first one is the ethical issue concerning the destruction of human embryos to generate hESCs, although these embryos were doomed to destruction anyway 2 . The second barrier to their potential medical use is the need of Human Leukocyte Antigen (HLA)-compatible [ ] hESC lines. Both issues are circumvented by induced pluripotent stem cells (iPSCs) as they are directly reprogrammed from adult cells 3 , [ ] thus opening the way to the derivation of patient-specific pluripotent stem cell (PCS) lines 4 . Moreover, by carefully selecting donor [ ] cells, for instance from donors who are homozygous for HLA, it will be possible to set up PSC banks that cover a large part of the genetic diversity of a population 5 . As iPSC derivation allows the production of cells from potentially any tissue in an unlimited supply [ ] in vitro and from any genetic background, many applications can therefore be foreseen, particularly: (i) modeling of human development, in vitro (ii) modeling of human genetic diseases, (iii) supply of normal or diseased cells for drug testing, (iv) supply of cells for cell in vitro therapies, including patient-specific cells, and (v) rejuvenation of old/senescent cells for regenerative medicine 6 . In addition, the iPS [ ] technology can be mastered with a reasonable amount of training and costs and is now successfully used in many laboratories worldwide (for a discussion on the iPS technology and its applications we refer the readers to other recent reviews 3,7 ). However, several reports are [ ] raising concerns about the potential weaknesses of iPSCs, such as the appearance of genetic or epigenetic abnormalities linked to the process of cell reprogramming. By altering the iPSC biological behavior or threatening their therapeutic safety, these defects could undermine their use for modeling diseases and medical applications. This problem might tip the balance back towards ESCs. In this review, we will focus on the genetic alterations observed in iPSCs and ESCs and discuss the fact that cell reprogramming is only one of the possible sources of genetic variation in PSCs.

TRACKING THE GENETIC DAMAGE PRODUCED BY CELL REPROGRAMMING

The generation of iPS by cell reprogramming opens the way to two potential sources of mutations. First, insertional mutagenesis due to the use of integrative vectors for cell reprogramming is an obvious cause of DNA damage (Figure 1). In the initial reports, Southern blot analyses showed that each iPSC clone was characterized by the presence of more than 20 integration events (about half a dozen integration events for each transcription factor) 8 . Although aberrant transcription of an oncogene following the integration of a retroviral vector in [ ] its proximity has been already described in gene therapy trials 9 , similar examples of insertional mutagenesis with functional [ ] consequences on endogenous genes are not a typical finding in iPSCs. Nonetheless, the definitive insertion of reprogramming vectors in the iPSC genome is a serious concern and the reactivation of the c-MYC transgene might have fatal consequences, as illustrated by the development of tumors in germline-competent adult chimeras, which were obtained using iPSCs generated with retroviral vectors, and in their progeny 10 . Safer reprogramming strategies have been developed and are now available, such as the use of non-integrating viral [ ] vectors, or the direct delivery of reprogramming transcription factor proteins or synthetic mRNAs (discussed below).

Second, the reprogramming process itself could also cause genomic alterations. expense of DNA integrity 26 28 .

[ -]

NOT ALL iPS CELLS HARBOR GENETIC DAMAGE

These findings rise the possibility that cell reprogramming might cause genomic damage and that, as a consequence, all iPSCs harbor genetic defects that may compromise their use for disease modeling or cell therapy. Several studies challenge this view. The karyotype analysis of more than 1700 human ESC and iPSC cultures (40 ESC and 219 iPSC lines) from 97 investigators in 29 laboratories did not reveal major differences between ESCs and iPSCs 14 . The overall incidence of karyotype abnormalities was 12.5 in human iPSCs and [ ] %

12.9 and hESCs, suggesting that cell reprogramming is not a significant source of karyotype abnormality. On the other hand, there were % differences concerning the occurrence of some abnormalities, such as higher incidence of trisomy 12 in hESCs and of trisomy 8 in iPCs, and most strikingly almost 10 of trisomy 17 in hESCs, but none in iPSCs. Interestingly, in this very large data set, recurrent % translocations were not found, differently from what is observed in many cancers. This is, however, a retrospective analysis and it is not clear whether the observed differences relate to the methods used to generate these PSC lines or whether they only reflect biases in the in cell expansion protocols. vitro

The idea that subkaryotypic abnormalities are unavoidable during cell reprogramming was also recently questioned. Qinlan et al. used whole-genome paired-end DNA sequencing to compare three iPSC lines that were generated by reprogramming mouse embryonic fibroblasts (MEF) with four lentiviral constructs coding for OCT4, SOX2, CMYC and KLF4 under the control of a doxycycline inducible promoter 29 . These low-passage iPSC lines were previously used to generate viable mice in tetraploid embryo complementation assays, [ ] undisputedly demonstrating that complete reprogramming had occurred. Unexpectedly, the authors found only one or two de novo structural genomic variations in each iPSC line. As a matter of fact, the most important source of genomic alteration was the de novo insertion of endogenous retroviral elements from the mouse leukemia virus (MLV) family that most likely originated from the CF-1 MEFs which were used as feeder cells. In another study, Young et al. analyzed by whole genome sequencing iPSC clones that were obtained using an OCT4, SOX2 and KLF4 polycistronic lentiviral vector and their parental fibroblasts. Their findings suggest that the single-nucleotide variations and the rare structural variations (SV) detected in the iPSC lines reflected the parental genomic background rather than the consequences of reprogramming 30 . These findings have been extended in human iPSC generated by using the canonical captures pre-existing mutations in the parental cell line during cloning, which is a necessary consequence of reprogramming. These data " obtained by whole genome sequencing of mouse iPSC lines were extended to human iPSC lines that were obtained using non-integrating Epstein-Barr nuclear antigen-1 (oriP/EBNA1) episomal vectors. Although more than 1,000 single nucleotide variations (SNVs) were found in each iPSC line, very few were within coding sequences 32 . Among those, 50 were synonymous changes, whereas SNVs that [ ] % modify the coding sequence were not enriched for known cancer-associated genes. Moreover, no CNV was found in the analyzed iPSC lines.

Overall, the studies suggesting that cell reprogramming may be inherently linked to DNA damage are counterbalanced by other works

showing that the mutation status of iPSC lines is similar to that of their parental cells. Nevertheless the reports about the occurrence of genome alterations during cell reprogramming are alarming as they highlight DNA damage. These contradictory reports could be bona fide explained mainly by two reasons: (1) different reprogramming techniques and/or primary cell types could result in dramatic differences in genome integrity; (2) different cell culture conditions could differentially affect PSC integrity and therefore be an important in vitro confounding variable. For instance, comparison of sequencing data of different iPSC lines derived by different groups, as done by Quinlan et al, demonstrated the huge variability in mutation frequency of such cells, strongly suggesting that some reprogramming protocols/culture conditions might more easily result in DNA damage 29 . From the data published to date, it is not possible to conclude [ ] whether the source cell type, the nature of the reprogramming vector (for instance with viral integration versus without integration) or the choice of the reprogramming gene combination, are correlated with the incidence of genomic abnormalities. In one report, some of these factors could be evaluated but none of these factors influenced the results 20 but these results are too preliminary to reach any [ ] conclusion.

CELL CULTURE-INDUCED GENOMIC ALTERATIONS

It will not be possible to precisely quantify the genetic abnormalities generated by cell reprogramming without knowing the role played by the cell culture conditions in the mutational process. Indeed, a cell culture protocol that favors the appearance of DNA alterations in PSCs will also do so during cell reprogramming and during the first passages which are required to amplify newly derived human iPSCs. In this regard, much can be learned from ESCs. ESCs are karyotypically normal at derivation 33,34 and the occurrence of sub-karyotypic abnormalities even faster when using single-cell dissociation passaging (unpublished data). These findings in ESCs indicate that the culture conditions used for cell reprogramming and cell expansion are a critical factor for DNA integrity during iPSC generation. Moreover, as several studies on genomic alterations in iPSCs were based on DNA samples isolated from cells before passage 20, it is not possible to exclude a contribution of the cell culture conditions to the DNA abnormalities that were ascribed to cell reprogramming 22 . Many issues need to be addressed to prevent these deleterious mutational events, particularly: (1) to identify the [ ] optimal culture conditions for PSC expansion and (2) to understand the mechanisms that drive the generation of genomic alterations.

Therefore, a quantitative assay is required to measure the genomic integrity of PSCs in order to compare the impact of different parameters on the incidence of mutational events.

IMPROVING CELL REPROGRAMMING

Overall, the genome integrity of human iPSCs and ESCs is comparable, although some differences have been reported. We will focus here on the measures that can be taken to diminish DNA damage during cell reprogramming.

To avoid insertional mutagenesis, several non-integrating approaches are currently available (see 3 for a discussion on these [ ]

The other concern is reprogramming-induced DNA damage. As its causes are still not clearly understood, the solution of this type of damage is still speculative. One direction could be to improve the reprograming transcription factor cocktail. Combinations that may improve cell reprogramming could associate the classical reprogramming factors with genes coding for components of the ATP-dependent BAF chromatin-remodeling complex 51 , or the use of the PSC-specific microRNAs miR-302 and miR-372 52 . To reduce PSC DNA integrity is still mainly assessed by karyotype analysis. Other approaches have been tested to overcome the obvious resolution limitations of the classic karyotyping techniques, for example CGH arrays or SNP microarrays; however, there is no consensus on the method to use to discriminate between the really worrying, possibly carcinogenic mutations and the DNA modifications with no or barely any impact on the biological behavior of PSCs, or the simple polymorphisms. For instance, in the preliminary report on the only published clinical trial based on human PSCs, the DNA integrity of the ESC line was confirmed by karyotyping of 20 G-banded metaphase cells and FISH analysis of chromosomes 12 and 17 in 200 interphase nuclei 54 .

[ ]

As highlighted by Prockop and Keating, the assessment of genomic stability in cultured cells is an old issue that was raised more than 50 years ago for cells in culture 55,56 . Indeed, genomic instability is a general feature of cultured cells and has been observed in many of the whole PSC genome will be available at a reasonable price.

PSCs without genome damage are obviously required for tissue modeling or medical applications. Next-generation sequencing technologies may be the solution to demonstrate that PSC lines are free of deleterious mutations 25 . However, a 100 undamaged DNA [ ] %

may not be necessary. The notion that only cells with pristine genetic information should be considered for use is certainly a preconception. Not only because the massive occurrence of polymorphisms in normal cells precludes a strict definition of pristine genetic information, but also because of the cell/tissue DNA information heterogeneity. DNA modifications arise continuously from early embryo development to adulthood and their occurrence increases during aging, leading to various degree of cell mosaicism 59 61 . The causes of [ -] these lesions are multiple and include, but are not limited to, random events promoted notably by reactive oxygen species (ROS), retrotransposition, UV (for skin fibroblasts) or senescence. Cell mosaicism can be readily observed in primary fibroblasts used for cell reprogramming when these cells are analyzed by using high-throughput techniques 22, 31 . As this genetic heterogeneity is uncovered by

[ ]
whole genome sequencing, there is an acute need to set up a demarcation line between acceptable and off-limit DNA alterations in the primary parental cells as well as in the freshly reprogrammed PSCs and long-term cultured PSCs. To this aim, the lesions that alter the cellular physiology or that pose a risk of malignant transformation must be identified 19 . As there is no means to assess the functionality Mutations that target a coding region and that are predicted to alter the protein sequence by giving rise to a new stop codon, a frameshift or a new splicing event.

Genetic defects that occur recurrently in cultured PSCs. To this aim, the creation of a PSC sequence database should be a key objective of scientists working on human PSCs. These recurrent mutations could be driving mutations (as opposed to passenger " " " " mutations) that confer a growth advantage to the cells. In addition, integration of these data with metabolic and phenotypic information will be instrumental for hierarchizing the abnormalities. Validation of these recurrent genetic abnormalities on large series of samples by independent techniques, for instance by digital droplet PCR for CNV, will be necessary 31 . In addition, measuring the impact of these

[ ]
abnormalities on the expression of the genes present in the vicinity by expression analysis will help in determining their functional consequences on the biology of PSC.

High-resolution, whole-genome sequencing data are ideally needed for the implementation of these guidelines, but they are not routinely available at the moment due to the actual costs of sequencing and the still complex pipeline analysis of such massive amount of data. Meanwhile, karyotyping will still be the reference test together with FISH analysis of chromosome 12, 17 and 20, and/or low-resolution CGH to exclude gross karyotype abnormalities.

CONCLUSION

The genetic integrity of iPSCs has been questioned and invoked as an argument against their medical use. However, these results have been challenged by reports that described the generation of iPSCs with low incidence of DNA sequence variations. Two factors may be central to PSC genetic stability: (i) the culture conditions, which can rapidly affect PSC genomic stability and (ii) the reprogramming protocols. As specific culture conditions are associated with very different rates of genetic variations, it can be anticipated that different starting cell types, reprogramming transcription factors and vectors, reprogramming protocols and early culture/passaging conditions can specifically affect the type/load of genomic damage linked to cell reprogramming. We therefore call for a careful quantification of the impact on genome integrity of all these factors in order to optimize the reprogramming protocols and ultimately provide safe iPSCs, both for research and clinical applications.

  [ ] retroviral method. The authors analyzed the iPSC lines by whole genome sequencing and digital droplet PCR and found that most of the iPSC CNVs reflect somatic mosaicism in the human skin 31 . These observations are in line with the hypothesis that cell reprogramming [ ] "

  has been thoroughly investigated. Aneuploid hESC clones can appear during cell culture. Such chromosomal abnormalities are often recurrent and gains of chromosome 12 (most frequently 12p), 17 (particularly 17q), 20 (most frequently 20q) or X have been often detected35 37 . These abnormalities might bestow a growth advantage and, as a consequence, cells carrying mutations[ -] are often selected during cell culture. The biological repercussions comprise increased proliferation capacity and decreased growth factor-dependence as well as modifications of the physiological differentiation capacities of such cells 38 . It must be noted that trisomy of [ ] chromosome 12, 17 or X is also detected in germ cell tumors and embryonal carcinoma cells. The same survival and/or cell proliferation pathways that are involved in germ and pluripotent cell tumor development could thus be involved also in the selection of aneuploid hESC clones during cell culture. As reported above, the extensive study of 40 hESC lines in which 1163 karyotypes were analyzed concluded that 12.9 of the hESC cultures displayed chromosomal aberrations 14 . This frequency could still increase if more sensitive methods % [ ] were used, such as microarrays 19,39,40 . Several factors could influence the incidence of these genetic aberrations, particularly the [ ] duration and the type of culture conditions. Enzymatic passaging might favor clonal selection of aneuploid clones, but this hypothesis has not been formally proven 39,41 . We recently observed that karyotypic abnormalities in hESCs can occur rapidly within ten passages and [ ]

  They include transient transfection 42 or the use of adenoviral vectors 43 , non-integrating viral vectors such as the [ ] [ ] oriP/EBNA1 episomal vectors 12 , non-integrating lentiviral vectors 44 and temperature-sensitive Sendai virus vectors 45 . Another removing the reprogramming vectors in the selected iPSC clones by using the piggyBac 46 or the Cre-LoxP excision " " [ ] system 47 after a transient integration step. Nevertheless, even when the vectors do not usually integrate into the DNA, or are excised [ ] after cell reprogramming, their insertion in the genome cannot be totally excluded. Thus, techniques that avoid DNA should preferred, such as the direct delivery of reprogramming cell-permeant proteins 48 , or of synthetic modified mRNAs 49 , which seem to have higher [ ] [ ] reprogramming efficacy. Direct protein delivery is less efficient probably because the innate immunity program needs to be activated via induction of the TLR3 pathway for effective cell reprogramming 50 . Indeed, in the classical reprogramming protocols, TLR3 is activated [ ]by the viral particles. Interestingly, addition of polyinosinic-polycytidylic acid (Poly I:C) to the culture medium as an alternative non-viral approach to activate TLR3, significantly improves and accelerates cell reprogramming of both mouse and human cells by cell-permeant proteins 50 .

  replication stress, c-MYC should be avoided. Of note, strategies that employ genes coding for viral oncoproteins (SV40) or inhibition of the p53 or p16 tumor suppressor pathways are clearly inacceptable when trying to minimize DNA damage 25 . As [ ] somatic-cell nuclear transfer and cell fusion, two other cell reprogramming techniques, are considered more efficient than the iPS technique 53 , these approaches might shed light on the molecular mechanisms underlying cell reprogramming and serve as a blue print to [ ] improve iPSC generation protocols. Finally, culture conditions are critical for the maintenance of PSC genome integrity. It will therefore be important to focus on how to improve PSC expansion . in vitro SETTING GUIDELINES TO ASSESS iPS CELLS Currently, there is no consensus on the extent of DNA damage that can occur during PSC generation/culture, the underlying mechanisms and the methods to avoid it. Therefore, a technique to analyze DNA integrity in PSCs and guidelines to interpret DNA variations are mandatory.

  including adult stem cells 57 . Yet, standards or guidelines about what should be considered as acceptable genome integrity of [ ] normal cells in culture are still lacking. As DNA sequencing technologies and their resolution (whole genome maps at single-base resolution) are improving very fast and their price rapidly decreasing 58 , it is time to anticipate the near future when the routine analysis [ ]

Figure 1

 1 Figure 1 Causes of genome alterations in iPSCs Genome alterations in iPSCs can be the consequence of: (A) variations or mutations that are present at high or low frequency in the parental cells; clonal selection during iPSC generation can select and/or fix these variants; (B) insertional mutagenesis, replicative stress, caused for example by forced expression of an oncogene, can result in DNA damage that might not be properly repaired; (C) iPSCs could be intrinsically genetically unstable, and sub-optimal culture conditions might favor the generation and/or the selection of clones that are fitter for in vitro culture conditions.

  Although iPSC lines have often normal karyotypes at early passages 11,12 , recent reports show that iPSCs, like ESCs, can display chromosomal aberrations 13 16 . The proportion of

	[	]	[ -]
	karyotype abnormalities may vary between iPSCs and ESCs, but most aberrations can be found in both PSC types 14 . In addition, [ ]
	detailed analyses, by using CGH microarrays 17,18 , SNP microarrays 19,20 or next-generation sequencing techniques 21,22 , suggest [ ] [ ] [ ]
	that more subtle abnormalities, such as copy number variations (CNV) and mutations, occur in iPSCs at much higher frequency than
	originally thought. These findings have raised serious concerns about the functional reliability and safety of iPSCs 23,24 . The involved [ ]
	mechanisms are multiple. The reprogramming process, particularly the forced expression of oncogenes, such as c-MYC which is often
	included in the reprogramming gene cocktail 18 , could cause replicative stress and consequently DNA damage. This can be repaired by [ ]
	the pathways which are normally used in metazoan cells to avoid malignant transformation, but if the repair machinery is defective,
	genetic aberrations might arise 25 . Indeed, inhibition of the p53 or p16 pathways greatly improves reprogramming efficiency, but at the [ ]

  of every known gene or to test the tumorigenicity of cultured cells, the only way is to make some assumptions about what kind of abnormalities could affect the physiology or safety of PSCs. Ideally, a PSC line should not display:

	Karyotype abnormalities, either in number or structure
	Genetic defects that have been reported in cancer cells, such as activating mutations of known oncogenes or inactivation of tumor
	suppressor genes (the COSMIC catalogue of somatic mutations in cancer database may help) 62 [ ]

[ ]

Acknowledgements:

We thank E. Andermarcher for critical reading of this manuscript. Supported by grants from the Association Fran aise contre les Myopathies

(ANR-08-BIOT-012-02).

Footnotes: DISCLOSURE OF POTENTIAL CONFLICTS OF INTEREST

The authors indicate no potential conflicts of interest.

LIST OF ABBREVIATIONS