Mohamed Mourad Hafidhi
email: mohamed.hafidhi@univ-ubs.fr

Emmanuel Boutillon
email: emmanuel.boutillon@univ-ubs.fr

Chris Winstead
email: chris.winstead@usu.edu

Reducing the impact of internal upsets inside the correlation process in GPS Receivers

This paper 1 examines two approaches to deal with internal logic upsets inside correlation process used in the tracking process of GPS receivers. These upsets can be produced due to process/voltage and temperature variations coupled with increased advancement of CMOS technology. If any upset occurs when computing the correlation function during each 10 ms, then errors are propagated in tracking loops, resulting in a loss of the GPS signal tracking and a distorted position given by the receiver. Results of experiments using a GPS receiver design are presented in this paper to evaluate the performance of each method. The two proposed solutions (the Feedback freezing loop (FFL) and the Last Correct Value (LCV) methods) offer a big interest compared to the classical Triple Modular Redundancy (TMR) method since they provide the same performance as the TMR with low area complexity. This work can be extended to any system using feedback loops information.

I. INTRODUCTION

Efforts are actually made to minimize power consumption and maximize battery life in mobile devices. In contrast operating at minimal power levels can be a source of internal logic upsets resulting in momentary/persistent faults in the device's behavior. This problem increases with the advancement of CMOS technology [START_REF] Jayanth Srinivasan | The impact of technology scaling on lifetime reliability[END_REF] in combination with process/temperature and voltage (PVT) variations [START_REF] Borkar | Tackling variability and reliability challenges[END_REF].

Global Positioning System (GPS) receivers are heavily used in mobile contexts, and there is motivation to minimize power consumption in these devices. Fault tolerance has been studied for data fusion in navigational systems that incorporate GPS receivers along with other sensor devices. For example, the authors of [START_REF] Jaradat | Enhanced, delay dependent, intelligent fusion for ins/gps navigation system[END_REF] consider an intelligent data fusion system to compensate for delayed or erroneous GPS data due to acquisition loss or signal obstruction. In [START_REF] Mohamed Mourad | Reliable gold code generators for gps receivers[END_REF] strategies have been proposed and evaluated for upset protection in Gold sequence generators in GPS receiver. [START_REF] Mourad Dridi | Reliable nco carrier generators for gps receivers[END_REF] presents four approaches for upset protection in the Numerically controlled oscillator (NCO) carrier generator for GPS receiver. To track satellites' signals, GPS receivers has to determine over time the carrier frequency and the distance satellite-receiver [START_REF] Tsui | Fundamentals of Global Positioning System Receivers[END_REF]. This is accomplished by comparing a locally generated signals 1 This work has received a French government support granted to the COMIN Labs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference ANR-10-LABX-07-01. It has also received support from the Brittany Region. Support was also from the US NSF award ECCS-0954747 and the Franco-American Fulbright Commission.

(Coarse/Acquisition codes and carrier) against the received signal, then moving smoothly in time local signals until it is time-aligned with the received signal. A correlation function is computed every 10 ms period to compare local signals with incoming signals. A maximum correlation output is achieved when the two signals are time aligned. Upsets when computing the correlation function can produce errors at its output. Because of feedback loops these errors will propagate over time and will corrupt the time estimation. Loss in the signal tracking process can be reported, forcing the receiver to restart the initial signal acquisition procedure. So it is increasingly important to deal with the impact of the upsets when they are generated in the correlation process.

The remainder of this paper is organized as follows. Sec. II gives detailed information in a compact form about GPS receiver architecture and rules of integrator modules. Sec. III describes some mechanisms to detect occurrences of upsets inside integrator modules. and Sec. IV presents the two proposed methods to reduce upsets' impact. Sec. V provides performance results, compares the probability for a receiver to be far from the reel (reference) position. Sec. VI offers conclusions.

II. GPS SIGNAL PROCESSING

A GPS is a well known technology that allows determining both the physical position and the absolute time of a receiver. The position in time and in space is determined thanks to a precise distance measurement with at least four GPS satellites. Each GPS satellite transmits a navigation message at 50-bits/s using the CDMA (Code Division Multiple Access) technology. The analytical expression of the transmitted signal of a satellite a is:

e a (t) = c a (t) d a (t) e 2πjf L1 t
where:

• d a (t): navigation message of the a th satellite,

• c a (t): a th Coarse/Acquisition (C/A) satellite code with a Binary Phase Shift Keying (BPSK) modulation (i.e. c a ∈ {-1, 1}), • f L1 : the carrier frequency in the L1 GPS Band (Open Sevice). GPS receiver has to demodulate the navigation message of different satellites in view to make the distance measurement. This involves two essential and sequential process: the acquisition process and tracking process. The acquisition process is the process by which the receiver identifies which satellites are in view. It is a a three-dimensional search to determine the GPS satellite identifier (which is the index of its associated C/A code), the code phase (represented by τ ), and the carrier frequency offset due to Doppler effect (represented by f d ).

Once estimated, the values of these parameters initiate the tracking modules. Since satellites are in continuous motion, the distance between any satellite and the receiver is dynamic. Besides to that, the carrier frequency of the received signal is also constantly changing in time due to Doppler shifts. Therefore, once acquired, GPS signals have to be tracked over time.

The GPS digital incoming signal is first demodulated by a frequency F L1 -F I (F I represent the intermediate frequency), and, after a pass band filter, sampled at a frequency of 4 MHz. The n th received sample s(n) is given by:

s(n) = k∈{1..K} A k (nT ) d k (nT -τ k (nT )) c k (nT -τ k (nT )) e jθ k n+ϕ k + w(n) (1)
where:

• θ k = 2πT f d(k) (nT ) + F I , • A k : Attenuation parameter of k th satellite signal,
• τ k : Propagation delay between satellite k and receiver,

• f d(k) : Doppler frequency of satellite k, • {1..K}: Set of satellites visible at time n,

• T : Sampling period (T = 0.25 µs),

• w(n): additive white gaussian noise.

As mentioned earlier the objectives of any tracking module is to extract the navigation message of the corresponding satellite. For each GPS satellite acquired, the GPS receiver associates a tracking channel module. A simplified representation of the channel tracking module is given in Fig. 1. Using estimated parameters fd and τ computed during the acquisition process, the GPS receiver produces, first, a local copy of the code C/A and the carrier frequency. Then a correlation process between locally generated signals and incoming signals is launched to compute the correlation function between signals during each 10 ms. C/A codes belongs to the family of Gold Pseudo Random Noise (PRN) sequences. Hence codes of all GPS satellites are orthogonal. Then, 10 ms-integrator output for the a th satellite is defined as, da (l)

10 ms s(n) c a (nT -τa (nT )) e -j( θan+ φa) (2)
Thus, a perfect time alignment (τ a = τ a ) and a perfect frequency and phase alignment ( θa = θ a and φa = ϕ a ) implies, according to (1) and ( 2), that da (l) = d a (l). The summation over 10 ms in (2) corresponds to the accumulation of 40,000 values (4 MHz clock during 10 ms) 2 . Based on the integrator output (which represents the result of the correlation function during 10 ms), discriminators and filters update the estimate Doppler frequency ( fd ) and the delay (τ ) in order to trade their evolution in time.

III. UPSETS DETECTION METHODS

Process/voltage/temperature (PVT) variations, coupled with continuous technology evolution and smaller transistor features have made internal logic upsets on devices more frequent. Logic upsets can manifest as momentary faults in the devices behavior, or as persistent faults. Occurrence of upsets can significantly affect the performance of circuits and can change output states. Therefore it is increasingly important to protect modern systems from upsets' impact to achieve acceptable reliability levels and maintain acceptable complexity and cost.

Detecting faults in combinatonational and sequential circuits is critical. Several techniques zre suitable to detect fault in combinatory circuit. First, the detection is done by the double modular redundancy (DMR) method. In this method the circuit is duplicated and a comparator is added to detect the incoherence between outputs of the two copies. In this case, area and power consumption increase 100 %. Besides to that DMR can detect upsets that affect only one copy.

To deal with these drawbacks, other design solutions for achieving transient-error detection were proposed. Detection can be realized by the double sampling method cited in [START_REF] Nicolaidis | Time redundancy based soft-error tolerance to rescue nanometer technologies[END_REF] and [START_REF]Circuit logique protege contre des perturbations transitoires[END_REF]. The main idea of this method is to take advantage of the temporal nature of transient faults, and achieve transient faults detection by using time redundancy. As shown in fig. 3, transient faults detection is realised by observing the output signals of a given circuit at two instants (double sampling).

The schema is composed of:

• A latch or flip-flop at the output of the circuit to make a delay version of the input signal. • A comparator to check the state of the original and the delayed versions. • An error detection flip-flop using as latching event the rising edge of the clock.

Simulations and performances are presented in [START_REF] Lorena | Cost reduction and evaluation of a temporary faults detecting technique[END_REF]. Results show that this new method uses a low hardware cost to guarantee a complete timing error and very soft errors detection. The upsets' detection can also be achieved by bulk built-in current sensors (BICS) [START_REF] Vargas | Seu-tolerant sram design based on current monitoring[END_REF], [START_REF] Balkaran | An efficient bics design for seus detection and correction in semiconductor memories[END_REF]. Based on the same approach (comparing the output at two different instant), BICS detect the ionisation current in the bulk by sensing the current. Two BICS are used in any solution to detect upsets: PBICS and NBICS. PBICS is used to detect output transitions from 0 to 1 and PBICS to detect output transitions from 1 to 0. Solutions based on the BICS approach has been proposed to make the detection in sequential and in combinational logic. The BICS appears in [START_REF] Daniel Gomez Toro | Soft error detection and correction technique for radiation hardening based on c-element and bics[END_REF] and [START_REF] Egas Henes | Using bulk built-in current sensors to detectsofterrors[END_REF] where a memory cell and 4-bit multiplier are respectively circuit under test. As for the double sampling method, with BICS approach we can also reduce area and power dissipation penalties compared to the DMR. To summarize, many are upsets' detection method that can be used to detect occurrence of upset in any application. In our case we will focus only on how we can limit upsets' impact after detecting their presence.

LIMIT UPSETS' IMPACT INSIDE CORRELATION PROCESS

In this section, we first define the method to assess the quality of a GPS with a noisy integrator. Then, assuming that error can be detected (see previous section), we define two methods to reduce efficiently the impact of errors.

A. Quality of a noisy GPS receiver

The methodology used to evaluate the quality of the proposed method is three fold:

• Use the front end of a GPS receiver to store incoming GPS signal s(n) over a significant period of time (few minutes). • Process the stored input GPS signal with a reference noiseless GPS receiver [START_REF] Dion | Personal communication[END_REF] to generate the set of successive estimated positions X(p) p=1...Q (The reference position). • Replace the noiseless GPS receiver with a noisy GPS receiver and process a Monte-Carlo simulation processing N times the stored incoming GPS signal to generate sets of noisy position Xi (p) p=1..Q , i = 1..N .

Then, we can study the link between the model of error injection, the hardware used to mitigate errors and the statistical difference between GPS noiseless position X and the GPS noisy position X.

1) Model of noise:

The model of noise used in this paper is very simple. For each active tracking module, every output of the integrator is assume to be exact with a probability (1 -p) or to be faulty with a probability p. In case of faulty result, a random value is uniformly drawn between the smallest and highest value of dk (l) to replace the exact value.

2) Metric to evaluate the quality of the noisy GPS receiver: Two metrics are used to evaluate the quality of the noisy GPS receiver. The first one is the mean squared error (EQM) defined by:

EQM (p) = E(| X -X | 2 ) ( 3
)
where p is the probability of the upsets' injection model (noise) at the output of the integrator 10 ms. This parameter gives us an idea about the distance, on average, we will be far from the reference position. The second parameter is:

f (d, p) = P r | X -X |> d (4)
Given the probability of error injection, p, f (d, p) reflect the probability that the noisy GPS receiver compute a position at a distance d from the reference position.

B. Feedback freezing loop (FFL) method

The first method we propose is to freeze the feedback loops once we detect occurrence of upset. The idea is that it is less harmful to inject nothing in the feedback loop than a wrong value. Fig. 2 illustrates the configuration used in this method. So if we detect that upsets has been produced during the 10 ms period (T 1 ), feedback loops are frozen (or blocked) until a new computation of the correlation function is achieved (i.e until the 10 next millisecond T 2 ). Code and carrier generators will use the same codes and carrier that have been generated during the previous 10 ms period (T 0 ). With this approach we can guarantee that errors in the correlation process don't propagate over time. Fig. 2 describes the configuration used for this method.

C. Last Correct Value (LCV) method

The second proposed method uses the last correct value of the integrator output to feed the feedback loop once we detect that an error has been produced when calculating the correlation at the next 10 ms-period T 1 . Fig. 2 describes the configuration used for this method. If da (l) is the result of the correlation function at time step l, at time step l + 1 and if errors have been produced when computing da (l + 1), the system uses the value da (l)to be sent to feedback loops instead of da (l + 1).

D. Triple modular redundancy (TMR)

TMR is a classical solution for fault tolerance in electronic systems [START_REF] Lyons | The use of triple-modular redundancy to improve computer reliability[END_REF]. It has been widely used for electronics in the space and nuclear industries to ensure reliable operation [START_REF] Cai | Implementation of a reconfigurable computing system for space applications[END_REF]. In a TMR system, the original module (the correlation function computaion (CorF) module) is replicated three times, and error correction is achieved by a majority vote operation. We note that there are several configurations, such as one-voter TMR, the triple-voter TMR or the restorative-feedback TMR [START_REF] Winstead | Error correction via restorative feedback in M-ary logic circuits[END_REF]but to minimize overhead we only consider the most basic TMR strategy as represented in fig. 4. If an upset occurs in any one of the three CorF modules, the other two CorF can correct and mask the fault assuming that the vote operation is not faulty. However once the voter is faulty, errors can be generated. These errors will propagate over time resulting in a fault at the output of the integrator 10 ms in Fig. 1.

V. PERFORMANCE RESULTS

All of the described methods were evaluated in term of EQM and Probability of an alignment of the reference position as defined in (3) and (4). Although they are simpl to implement and with a very low cost area, LCV and FFL methods offers a high gain performance. Fig. 5 compares EQM of different methods. In this figure, EQM = 4.65 10 6 m (which represent the reference position) is equivalent to a loss of the GPS signal tracking. As we can see in Fig. 5, for a probability of upsets up to p = 0.4, LCV method gives an error in the position that does not exceed on average 2 m. Thus the LCV methods permit to achieve an acceptable reliability even with high upset probability. Fig. 6 and Fig. 7 show the variation of this probability as a function of the distance between the noisy and the reference position when p = 0.1 and p = 0.2 respectively. For instance when p = 0.1, the probability of been 1 m far from the reference position goes from 0.6 to 0.009 using the LCV method.

VI. CONCLUSION

Our results show an improvement in the probability of the removal from the reference position and in term of EQM with the LCV and FFL methods. Up to a probability of upset p = 0.4, with LCV configuration we guarantee an acceptable reliability (EQM ¡ 2 m) for a low hardware complexity compared to the TMR strategy.

Fig. 1 .

 1 Fig. 1. Top Level of the simplified tracking channel module

Fig. 2 .

 2 Fig.2. Configuration of proposed method. We add to the original tracking channel either the FFL or the LCV configuration

Fig. 3 .

 3 Fig.3. Double sampling approach[START_REF]Circuit logique protege contre des perturbations transitoires[END_REF]

Fig. 4 .

 4 Fig. 4. TMR strategy of protectionIn the next section, simulation results show the quality of the two proposed methods compared to the classical TMR.

Fig. 5 .

 5 Fig.5. EQM = f (p); When p = 10 -1 , EQM goes from 1000 km to an orror in the position of about 1 meters (using FFL method) and few centimeters (with LCV method).

Fig. 6 .

 6 Fig.6. The probability of the removal (in meters) from the reference position when the upset probability p = 0.1

Each bit of the navigation message have a duration of 20 ms (50 Hz bit rate). Thus, assuming perfect bit timing synchronisation, da(l) is constant while computing[START_REF] Borkar | Tackling variability and reliability challenges[END_REF]

× Incoming signal × Integrator 10 ms Code generator Carrier NCO Code Filter Carrier Filter Code Discriminator Carrier Discriminator


Mourad Dridi
email: mounad.dridi@univ-ubs.fr

Mohamed Mourad Hafidhi #

Chris Winstead
email: chris.winstead@usu.edu

Emmanuel Boutillon
email: emmanuel.boutillon@univ-ubs.fr

Reliable NCO carrier Generators for GPS Receivers

Four strategies for upset protection in NCO (Numerically controlled oscillator) generators are presented and compared in this paper 1 . This work is motivated by the increased prevalence of single event upsets due to process/voltage/temperature variations coupled with the increased advancement of CMOS technology. The occurrence of these upsets in the GPS context (the system under the test in this paper) can lead to a loss of the GPS signal tracking, then the corresponding satellite signal must be re-acquired resulting in high energy expenditure and time delay. It can result also in a corrupted position given by the GPS receiver. Results of experiment using a GPS receiver design are presented in this paper to evaluate the performance of proposed methods. This work can be extended and generated to any system using feedback loops information.

I. INTRODUCTION

During the last decades, increased advancement of CMOS technology coupled with process/voltage/temperature (PVT) variations, have made internal upsets on devices more frequent [1], [START_REF] Borkar | Tackling variability and reliability challenges[END_REF]. These upsets manifest as momentary fault in the device's behavior and their occurrence can significantly affect the circuit performance and output states. Therefore, it is extremely important to protect systems from fault's impact to achieve acceptable reliability and maintain low complexity and cost. Fault tolerance has been studied for mobile systems to build mobile devices more robust to transient errors due to PVT variations. For example, the author of [START_REF] Gimmler-Dumont | A cross-layer reliability design methodology for efcient, dependable wireless receivers[END_REF] proposed some resilience techniques for double-iterative multiple-input multiple-output bit-interleaved coded modulation (MIMO BICM) receivers. These techniques have been realized on different layers of abstraction in order to have a mutual trade-off of system performance, hardware reliability and implementation complexity.

Global Positioning System (GPS) receivers are extensively used in mobile context, and there is motivation to minimize power consumption and maximize battery life in these devices. Fault tolerance has been studied for data fusion in navigational systems that incorporate GPS receivers along with other sensor devices. For example, the authors of [START_REF] Jaradat | Enhanced, delay dependent, intelligent fusion for ins/gps navigation system[END_REF] consider an intelligent data fusion system to compensate for delayed or erroneous GPS data due to acquisition loss or signal obstruction. In 1 This work has received a French government support granted to the COMIN Labs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference ANR-10-LABX-07-01. [START_REF] Mohamed Mourad | Reliable gold code generators for gps receivers[END_REF], strategies have been proposed and evaluated for upset protection in Gold sequence generators. Paper [START_REF]Reducing the impact of internal upsets inside the correlation process in gps receivers[END_REF] presents two approaches to deal with internal logic upsets inside correlation process used in the tracking process of GPS receivers.

GPS receivers use Numerically Controlled Oscillators (NCO) to produce a local copy of the carrier frequency offset (due to Doppler offset) of incoming satellite signals. When the NCO output is corrupted, distance measurements in the GPS receiver are disturbed and this leads to a faulty localization of the receiver. This problem can be mitigated if the NCO is protected. In this context, this paper presents four solutions to make the NCO module more robust to upsets and evaluates those solutions in terms of their performance. This work can be generated for any NCO in feedback loop systems. The remainder of this paper is organized as follows. Sec. II gives detailed information in a compact form about the GPS receiver architecture and rules of NCO generators. Sec. III presents the Triple Modular Redundancy (TMR) method. Sec. IV, V and VI describe three methods using the modular redundancy. Sec. VII provides synthesis and performance results. Finally Sec. VIII offers conclusions.

II. GPS SIGNAL PROCESSING

GPS is the US Global Navigation Satellite System (GNSS) that provides free positioning and timing service. Each GPS satellite broadcasts a navigation message at 50 bits/s using the Code Division multiple Access (CDMA) technology. The analytical expression of the transmitted signal of a satellite k is

e k (t) = c k (t) d k (t) e (2πjf L1 t)
, where:

• d k (t): navigation message of the k th satellite; • c k (t): k th Coarse/Acquisition (C/A) satellite code with a Binary Phase Shift Keying (BPSK) modulation (i.e. c k ∈ {-1, 1}); • f L1 : the carrier frequency in the L1 GPS Band (Open Sevice).
At least four GPS satellites are needed for a GPS receiver to be able to compute the position in time and space. The incoming digital GPS signal in the receiver is first demodulated by a frequency F L1 -F I (F I represents the intermediate

s(n) = k∈{1..K} A k (nT ) d k (nT -τ k (nT )) c k (nT -τ k (nT )) e jθ k n+ϕ k + w(n) (1)
where:

• θ k = 2πT f d(k) (nT ) + F I , • A k : Attenuation parameter of k th satellite signal,
• τ k : Propagation delay between satellite k and receiver,

• f d(k) : Doppler frequency of satellite k, • {1..K}: Set of satellites visible at time n,

• T : Sampling period (T = 0.25 µs),

• w(n): additive white Gaussian noise. For each GPS satellite visible in view, the GPS receiver tries to extract the corresponding navigation message. To do so, there are two essential and sequential processes: the acquisition and tracking process. The acquisition process is the process by which the receiver identifies which satellites are in view. It is a a three-dimensional search to determine the GPS satellite identifier (i.e the C/A code c k associated to the satellite), code phase (represented by τ ), and carrier frequency offset due to Doppler effect (represented by fd ). Since satellites are in continuous motion, distances between satellites and the receiver are dynamic. The carrier frequency of the received signal is also constantly changing in time due to Doppler shifts. Therefore, once the acquisition process is accomplished, GPS signals have to be tracked over time.

Each acquired satellite will be tracked by a unique tracking channel. A simplified representation of the channel tracking module is given in Fig. 1. The main objective of each channel tracking module is to extract the navigation message of the corresponding GPS satellite. To do that, using estimated parameters fd and τ , the GPS receiver produces a local copy of the C/A code (via the Gold Code generators) and the carrier frequency (via the NCO component). Then, a comparison between locally generated signals and incoming signals, using a correlation process and for a period of 10 ms, is made. C/A codes belongs to the family of Gold Pseudo Random Noise (PRN) sequences. Hence codes of all GPS satellites are orthogonal. Then, the 10 ms integrator output for the a th satellite is defined as, da ( )

10 ms s(n) c a (nT -τa (nT )) e -j( θan+ φa) (2)
Thus, a perfect time alignment (τ a = τ a ) and a perfect frequency and phase alignment ( θa = θ a and φa = ϕ a ) implies, according to (1) and ( 2), that da ( ) = d a ( ). The summation over 10 ms in (2) corresponds to the accumulation of 40,000 values (4 MHz clock during 10 ms) 2 . Taking in account the correlation function result, the estimate of the Doppler frequency ( fd ) and the delay (τ ) is made via discriminators and filters. The carrier loop filter and the carrier discriminator produce an accurate estimate of the Doppler offset of incoming signals while the code loop filter and the code discriminator produce an accurate estimate of the delay offset [START_REF] Tsui | Fundamentals of Global Positioning System Receivers[END_REF].

∀n ψ(n) = ψ (n -1) + θ ( n/N s ) ( 3
)
where x is the floor function of x and N s is total number of samples generated by the NCO in 10 ms (N s = 40000 in the GPS context). If the NCO outputs are disturbed due to logic internal upsets, then the corresponding carrier frequency becomes corrupted, resulting in a loss of the GPS signal tracking. This problem can be mitigated if the NCO is protected by an error corrector method. Next sections will present different methods that can be proposed to protect the NCO carrier from momentary errors.

III. TRIPLE MODULAR REDUNDANCY (TMR)

TMR is a classical solution for fault tolerance in electronic systems [START_REF] Lyons | The use of triple-modular redundancy to improve computer reliability[END_REF]. It has been widely used for electronics in the space and nuclear industries to ensure reliable operation [START_REF] Cai | Implementation of a reconfigurable computing system for space applications[END_REF]. In a TMR system, the original module (the adder + register in this case) is replicated three times, and error correction is achieved by a majority vote operation. We note that there are several configurations, such as one-voter TMR, the triple-voter TMR or the restorative-feedback TMR [START_REF] Winstead | Error correction via restorative feedback in M-ary logic circuits[END_REF]. In this paper we will use the one-voter and the triple-voter TMR as represented in Fig. 3. If an upset occurs in anyone of the three replicas, the other two replicas can correct and mask the error. But once two of the three replicas fail, an error will be generated in the system and propagated over time. We note also that errors in the voter are extremely harmful for the system because they will not be corrected. To deal with these drawbacks, the triple-voter TMR also triplicate the voters so that a vote error can be corrected in a subsequent circuit stage.

IV. RENDEZ-VOUS CHECKING METHOD (RVC)

In a TMR system, the original module is replicated three times and error correction is achieved by a majority vote operation. In this case area and power consumption increase of more than 300 %. To design a NCO tolerant to transient error with low complexity and cost we propose here to use the assumption that the NCO input remains constant during each 10 ms. So from (3), we can deduce that:

∀ ∈ N ψ(( + 1) N s ) = ψ( N s ) + N s θ( ) (4)
Based on (4), we add to the original architecture (i.e the adder and the register) a correction module composed by a multiplier and adder as configured in Fig. 4. At the end of each 10 ms, the correction module produces a new output (ψ c ) by substituting ψ c in place of ψ in (4). The Rendez-vous Checking Method will consist on a meeting (Rendez-vous) organised at the end of each 10 ms to check the coherence of values between ψ c and ψ o (where ψ o is the output of the original module given by substituting ψ o in place of ψ in (3) ). Once an incoherence (ψ o = ψ c ) is declared, it is ψ c that will be used later clocks.

Fig. 4. Rendez-vous checking method configuration Fig. 5 illustrates how we are able to correct errors occurring when the NCO is operating. However, when an error occurs when the correction module is operating, an uncorrectable failure results in the system. We have to note here that the probability to have a faulty correction module is low once the correction module output is computed during 10 ms. Fig. 5. Correction of errors in RvC method; when an error occurs at time step t1, it can be corrected at time at 20 ms thanks to the coorrect output of the correction module.

V. TIME FREEZING METHOD (TF)

In order to make the NCO carrier generator more robust against upset, the original module is duplicated. Then we use a comparator to detect incoherence between replicas output. Let's consider :

VIII. CONCLUSION

Our results show an improvement in the probability of the removal from the reference position and in term of EQM with each methods introduced previously. For a small upset probability, Rendez Vous Checking method provides an acceptable reliability for a low cost configuration. For a high upset probability, Robust Time Freezing and Time Freezing method provides the big interest for reducing errors in the NCO carrier.

[1] P. B. J. A. Jayanth Srinivasan, Sarita V. Adve, "The impact of technology scaling on lifetime reliability," The International Conference on

Fig. 1 .

 1 Fig. 1. Generic digital receiver channel block diagram

Fig. 2 .

 2 Fig. 2. NCO generator block diagramFig.2shows a block diagram of the carrier loop NCO and its sine and cosine mapping functions. The NCO carrier receives at its input the estimated Doppler offset value θ. The value of θ remains constant each 10 ms (tracking loops are operating at 100 Hz). If we consider ψ the output of the NCO carrier, ψ is produced each 0.25 µs (i.e at 4 MHz) and is defined as

  Fig. 3. NCO carrier with TMR protection

Fig. 8 .

 8 Fig.8. EQM = f (p); When p = 10 -1 , EQM goes from 50 meters to an orror in the position less than 1 meter using the TMR 3 Voter method.

Fig. 9 .

 9 Fig. 9. The probability of the removal (in meters) from the reference position when p = 10 -1

  No Correction p =0.01 TF method p =0.01 RTF method p =0.01 TMR 3 Voter p =0.01

Fig. 10 .

 10 Fig. 10. The probability of the removal (in meters) from the reference position when p = 10 -2

Each bit of the navigation message have a duration of 20 ms (50 Hz bit rate). Thus, assuming perfect bit timing synchronisation, da( ) is constant while computing[START_REF] Borkar | Tackling variability and reliability challenges[END_REF]

• ψ(n): the output of the NCO at time step n

Code Discriminator

Code Filter Code Generator

• ψ 1 (n): the output of Adder 1 at time step n • ψ 2 (n): the output of Adder 2 at time step n Once one of the two adders are faulty (i.e ψ 1 (n) = ψ 2 (n), e = 1), the correct previous output of the NCO (i.e ψ(n -1)) is used to feedback adders as shown in Fig. 6. Fig. 6. Time Freezing Method configuration However, if the register R is faulty, the TF configuration fails and errors will propagate over time in the tracking loops.

VI. ROBUST TIME FREEZING (RTF) METHOD

To deal with TF method drawbacks we propose in this method to triplicate the register as shown in Fig. 7. The register R in Fig. 7 is used to store the NCO output and to feedback adders with the last correct output value of the NCO when an incoherence between adders output is detected. As for the TMR 3 Voter, any simple error can be detected and corrected in the RTF configuration.

VII. PERFORMANCE RESULTS

In this section, we first define the method to assess the quality of a GPS with a noisy NCO carrier. Then, all of the described methods were evaluated and compared.

A. Quality of a noisy GPS receiver

The methodology used to evaluate the quality of the proposed method is threefold:

• Use the front end of a GPS receiver to store incoming GPS signal s(n) over a significant period of time (few minutes).

• Process the stored input GPS signal with a reference noiseless GPS receiver [START_REF] Dion | Personal communication[END_REF] to generate the set of successive estimated positions X(p) p=1...Q . • Replace the noiseless GPS receiver with a noisy GPS receiver and process a Monte-Carlo simulation processing N times the stored incoming GPS signal to generate sets of noisy position Xi (p) p=1..Q , i = 1..N . Then, we can study the link between the model of error injection, the hardware used to mitigate errors and the statistical difference between the noiseless GPS position X and the noisy GPS position X.

1) model of noise: The model of noise used in this paper is very simple. For each active tracking module, every output of the NCO module is assume to be exact with a probability (1 -p) or to be faulty with a probability p. In case of faulty result, a random random is uniformly drawn 0 and 2π to replace the exact value.

2) Metric to evaluate the quality of the noisy GPS receiver: Two metrics are used to evaluate the quality of the noisy GPS receiver. The first one is the Mean Squared Error (EQM) defined by:

where p is the probability of the upsets' injection model (noise) at the output of the integrator 10 ms. This parameter gives us an idea about how many meters, on average, the position given by the noisy GPS is far from the position given by the reference GPS. The second parameter is:

Given the probability of error injection, p, f (r, p) reflect the probability that the noisy GPS receiver compute a distance at a distance greater than r meter from the noiseless GPS receiver.

B. Evaluation Methods

All of the described methods were evaluated in term of EQM and probability of an alignment of the reference position as defined in ( 5) and [START_REF]Reducing the impact of internal upsets inside the correlation process in gps receivers[END_REF]. Fig. 8 compares the mean square squared error (EQM) of different method. As we can see from this figure (Fig. 8), Rendez-vous Checking Method provides us a very important result (EQM < 10 m) for the low error probability but it lose its efficiency for the higher error probability. Time Freezing Method and Robust Time Freezing Method correction is more efficient than the TMR one voter and they are very close to the TMR 3 voter. For example, for a probability of upsets up to p = 10 -2 , RTF method gives an error of the position of about 10 meters, 2 meters more than the EQM for the TMR-3 Voter method. In Fig. 9 and Fig. 10 we fix the probability of upsets (p = 10 -1 , for p = 10 -2 respectively) and we try to compute the probability to be far from the reference position. For instance when p = 0.01, the probability of been 10 m far from the reference position goes from 0.42 (No correction mechanism is added to the NCO) to 0.01 using the TF method.


Loïc Chaumont
email: loic.chaumont@univ-angers.fr

Thi Ngoc

Anh Nguyen
email: nguyen@math.univ-angers.fr

Thi Ngoc Anh

ON MUTATIONS IN THE BRANCHING MODEL FOR MULTITYPE POPULATIONS

Keywords: . 2010 Mathematics Subject Classification. 60J80 Multitype branching forest, mutation, forest of mutations, emergence time. LOÏC CHAUMONT

The forest of mutations associated to a multitype branching forest is obtained by merging together all vertices of its clusters and by preserving connections between them. We first show that the forest of mutations of any mulitype branching forest is itself a branching forest. Then we give its progeny distribution and describe some of its crucial properties in terms the initial progeny distribution. We also obtain the limiting behaviour of the number of mutations both when the total number of individuals tends to infinity and when the number of roots tends to infinity. The continuous time case is then investigated by considering multitype branching forests with edge lengths. When mutations are non reversible, we give a representation of their emergence times which allows us to describe the asymptotic behaviour of the latters, when the ratios of successive mutation rates tend to 0.

Introduction

The homogeneous multitype branching hypothesis provides a relevant model of population growth in the absence of any competitive or environmental constraint. In particular, it is widely used in population genetics, when studying successive mutations whose accumulation leads to the development of cancer. Then determining the statistics of the emergence times of mutations, or evaluating the distribution of the population size of mutant cells at any time become important challenges. In the extensive literature on the subject, let us simply cite [START_REF] Iwasa | Population genetics of tumor suppressor genes[END_REF], [START_REF] Haeno | The evolution of two mutations during clonal expansion[END_REF], [START_REF] Durrett | Evolution of resistance and progression to disease during clonal expansion of cancer[END_REF], and [START_REF] Durrett | Population genetics of neutral mutations in exponentially growing cancer cell populations[END_REF].

This work is concerned with the mathematical study of mutations in multitype branching frameworks. We first focus on the problem of the total number of mutations under very general assumptions. This number is not a functional of the associated branching process and its study requires the complete knowledge of the multitype branching structure, that is the underlying plane forest. Then we show that the forest of mutations associated to any multitype forest, is itself a multitype branching forest whose progeny distribution can be explicitely computed. This result allows us to investigate the asymptotic behaviour of the number of mutations, when either the total population or the initial number of individuals tend to infinity.

When time is continuous, we are mainly interested in emergence times of new mutations in the non reversible case. The characterisation of these times requires a good knowledge of the corresponding multitype branching process and the main tool in this study consists in a recent extention of the Lamperti representation in higher dimensions. Emergence times are then expressed in terms of the underlying multivariate compound Poisson process, which allows us to obtain some accurate approximations of their law.

We start with some preliminaries on the coding of multitype branching forests by multivariate random walks in Section 2.1. Then we state and prove our results on the A cluster or a subtree of type i ∈ [d] of a d-type forest (f , c f ) ∈ F d is a maximal connected subgraph of (f , c f ) whose all vertices are of type i. Formally, t is a cluster of type i of (f , c f ), if it is a connected subgraph whose all vertices are of type i and such that either the root of t has no parent or the type of its parent is different from i. Moreover, if the parent of a vertex v ∈ v(t) c belongs to v(t), then c f (v) = i. Clusters of type i in t 1 are ranked according to the order of their roots in the breadth first search order of t 1 , see Figures 1 and2. Then if the number of clusters of type i is finite in t 1 , we continue by ranking clusters of type i in t 2 , and so on. Note that with this procedure, it is possible that clusters of t k , t k+1 , . . . , for some k, are not ranked. We denote by t

k , . . . } is called the subforest of type i of (f , c f ). We denote by u (i) [START_REF] Alexander | Conditional distributions and waiting times in Multitype branching processes[END_REF] , u (i) 2 , . . . the elements of v(f (i) ), ranked in the breadth first search order of f (i) . The subforests of the 2-type forest given in Figure 1 are represented in Figure 2.

(2)

11 [START_REF] Iwasa | Population genetics of tumor suppressor genes[END_REF] 13 15 14 t 2

(2) t 3

(2) Figure 2. The subforests of the 2-type forest given in Figure 1 with their deapth first search labeling.

t 5 (1)
To any forest (f , c f ) ∈ F d , we associate the forest of mutations, denoted by ( f , cf ) ∈ F d , which is the forest of F d obtained by aggregating all the vertices of each subtree of (f , c f ) with a given type, in a single vertex with the same type, and preserving an edge between each pair of connected subtrees. An example is given in Figure 1.

For a forest (f , c f ) ∈ F d and u ∈ v(f ), when no confusion is possible, we denote by p i (u) the number of children of type i of u. For each i ∈ [d], let n i ∈ Z + ∪ {∞} be the number of vertices in the subforest f (i) of (f , c f ). Then let us define the d-dimensional chain x (i) = (x i,1 , . . . , x i,d ), with length n i and whose values belong to the set Z d , by x

(i) 0 = 0 and if n i ≥ 1, (2.1) x i,j n+1 -x i,j n = p j (u (i) n+1 ) , if i = j and x i,i n+1 -x i,i n = p i (u (i) n+1 ) -1 , 0 ≤ n ≤ n i -1 , where (u (i) n ) n≥1
is the labeling of the subforest f (i) in its own breadth first search order. Note that the chains (x i,j n ), for i = j are nondecreasing whereas

(x i,i n ) is a downward skip free chain, i.e. x i,i n+1 -x i,i n ≥ -1, for 0 ≤ n ≤ n i -1. Besides, if n i is finite, then n i = min{n : x i,i n = min 0≤k≤n i x i,i k }.
Let us also mention that from Theorem 2.7 of [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF], when trees of (f , c f ) are finite, the data of the chains x (1) , . . . , x (d) together with the sequence of ranked roots of (f , c f ), allow us to reconstruct this forest.

Let us now apply this coding to multitype branching forests. Let ν := (ν 1 , . . . , ν d ), where ν i is some distribution on Z d + . We consider a branching process with progeny distribution ν, that is a population of individuals which reproduce independently of each other at each generation. Individuals of type i give birth to n j children of type j ∈ [d] with probability ν i (n 1 , . . . , n d ). For i, j ∈ [d], we denote by m ij the mean number of children of type j, given by an individual of type i, i.e.

m ij = (n 1 ,...,n d )∈Z d + n j ν i (n 1 , . . . , n d ) .
We say that ν is non singular if there is i ∈

[d] such that ν i (n : n 1 + • • • + n d = 1) < 1. The matrix M = (m ij ) is said to be irreducible if for all i, j, m ij < ∞ and there exists n ≥ 1 such that m (n) ij > 0, where m (n)
ij is the ij entry of the matrix M n . If moreover the power n does not depend on (i, j), then M is said to be primitive. In the latter case, according to Perron-Frobenius theory, the spectral radius ρ of M is the unique eigenvalue which is positive, simple and with maximal modulus. If ρ ≤ 1, then the population will become extinct almost surely, whereas if ρ > 1, then with positive probability, the population will never become extinct. We say that ν is subcritical if ρ < 1, critical if ρ = 1 and supercritical if ρ > 1. We sometimes say that µ is irreducible, primitive, (sub)critical or supercritical, when this is the case for M . By multitype branching forest with progeny distribution ν, we mean a sequence with a finite (deterministic) or infinite number of independent multitype branching trees with progeny distribution ν. A multitype branching forest will be considered as a random variable defined on the probability space (Ω, F, P ) and with values in F d . To any multitype branching forest F , we associate the random sequences

X = {X (i) , i ∈ [d]}, where X (i) = {(X i,1 n , . . . , X i,d n ), 0 ≤ n ≤ n i },
which are constructed as in (2.1). It has been proved in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF], Theorem 3.1 that if F is a primitive and (sub)critical branching forest with a finite number of trees, then X (i) , i ∈ [d] are independent random walks whose step distribution νi is defined by

(2.2) νi (k 1 , . . . , k d ) := ν i (k 1 , . . . , k i-1 , k i + 1, k i+1 , . . . , k d ) , for all (k 1 , . . . , k d ) ∈ Z d
+ , and stopped at the smallest solution (N 1 , . . . , N d ) of the system

(2.3) x j + d i=1 X i,j (N i ) = 0 , j ∈ [d] .
In this equation, N i is the total number of vertices of type i in F and x i is the total number of trees in this forest whose root is of type i. We will say that F is issued from x = (x 1 , . . . , x d ). Note that the variables N i are random, whereas the x i 's are deterministic.

2.2. The total number of mutations and its asymptotics. A mutation of type i, is the birth event of an individual of type i from an individual of any type j = i. The aim of this section is to study the evolution of mutations in a multitype branching forest. Our main result asserts that the forest of mutations, that is the forest obtained by merging together all the vertices of a same cluster, is itself a branching forest if and only if for each i ∈ [d], one of the following conditions is satisfied,

(A i ) m ii ≤ 1 , (B i )
m ii > 1 and for all j = i, m ij = 0.

Moreover, its progeny distribution can be expressed in terms of this of the initial forest. Note that the branching property of the forest of mutations is intuitively clear. In the neutral case, it has been pointed out in [START_REF] Taïb | Branching processes and neutral evolution[END_REF].

Theorem

k i = 0}
, which is defined by

(2.4) µ i (k) = n≥1 n -1 ν * n i (k + (n -1)e i ) , k ∈ S i , if (A i ) is satisfied. If (B i
) is satisfied, then µ i is the Dirac mass at 0. Moreover µ satisfies the following properties:

1. Let M = ( mij ) be the mean matrix of µ and let r ≥ 1. Then µ i admits moments of order r if and only if either for all j = i, m ij = 0 or ν i admits moments of order r and m ii < 1. In the latter case, for all i, j such that i = j, mij =

m ij 1-m ii . 2. Assume that mij < ∞, for all i, j ∈ [d]. Then M is irreducible if and only if M is irreducible. If M is primitive, then so is M . The converse is not true. 3.
Assume that M is primitive, then M is subcritical (resp. critical, resp. supercritical) if and only if M is subcritical (resp. critical, resp. supercritical). If for some i ∈ [d], none of the conditions (A i ) and (B i ) holds, then there is j = i such that individuals of type i in F give birth to an infinite number of children of type j with positive probability. Therefore F is not a branching forest in our sense.

Proof.

Since the result only bears on the progeny law of forests, we do not loose any generality by assuming that F has an infinite number of trees. Then the stochastic processes X = {X (i) , i ∈ [d]} obtained from F , as in (2.1) are defined on the whole integer line {0, 1, . . . }. Note that their definition slighly extends the definition which is given in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF]. Indeed, without any more assumption on ν, trees of the forest can be infinite, so that the process X is not necessarily a coding of the forest, that is, if some trees are infinite then it is not possible to reconstruct the whole forest from X and the sequence of its roots. However, it is straightforward to check that X (i) , i ∈ [d] are independent random walks and that the step distribution of X (i) is νi , which is defined in (2.2). In particular, the law of X characterizes this of F . Now, let us consider the forest of mutations F . By construction, this forest is composed of an infinite number of independent and identically distributed trees. Hence, in order to show that F is a branching forest, it suffices to show that its trees are branching trees.

Let us denote by {X (i) , i ∈ [d]} the process which is defined from F as in (2.1). Let i ∈ [d] and assume first that (A i ) holds. Then we define the first passage time process of the random walks

X i,i , i ∈ [d] by, τ (i) k = inf{n ≥ 0 : X i,i n = -k} , k ≥ 0 .
Since m ii ≤ 1, then from the law of large numbers, lim inf n→∞ X i,i n = -∞, a.s., so that τ (i) k is almost surely finite for all k ≥ 0 and lim k→∞ τ

(i) k = ∞, a.s. Moreover, for all i, j ∈ [d], X i,j k = X i,j (τ (i) k ) , k ≥ 0 .
Indeed, the effect of the time change by τ (i) k is to merge all vertices of a same cluster of type i into a single vertex. Note that X (i) , i ∈ [d] are independent random walks.

Assume with no loss of generallity that the root of the first tree in F has type 1, then a slight extention Theorems 2.7 and 3.1 in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF] to any progeny distribution, allows us to show that this first tree is coded by the processes (X

(i) k , 0 ≤ k ≤ N i ), i ∈ [d], where (N 1 , . . . , N d ) is the smallest solution of the system (2.5) r j + d i=1 X i,j (N i ) = 0 , j ∈ [d] ,
and (r 1 , . . . , r d ) = (1, 0, . . . , 0). Note that in our case, N i can be infinite. This extended notion of smallest solution is defined in [START_REF] Chaumont | Breadth first search coding of multitype forests with application to Lamperti representation[END_REF], see Lemma 1 therein. This coding result implies that the first tree in F can be reconstructed from the processes (X

(i) k , 0 ≤ k ≤ N i ), i ∈ [d]
and applying part 3. of Theorem 3.1 in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF], we obtain that this tree is a branching tree whose progeny distribution µ = (µ i , i ∈ [d]) is given by

µ i (k 1 , . . . , k d ) = P (X (i) 1 = (k 1 , . . . , k i-1 , -1, k i+1 , k d )) , (k 1 , . . . , k d ) ∈ S i .
Then in order to make this law explicit in terms of ν, we apply the Ballot theorem for cyclically exangeable sequences due to Takács [START_REF] Takács | The probability law of the busy period for two types of queuing processes[END_REF]. Since conditionally on X i,j , i = j, X i,i is downward skip free with cyclical exchangeable increments, we have for all (k 1 , . . . , k d ) ∈ S i , P (X

(i) 1 = (k 1 , . . . , k i-1 , -1, k i+1 , . . . , k d )) = n≥1 P (X (i) n = (k 1 , . . . , k i-1 , -1, k i+1 , . . . , k d ), τ (i) 1 = n) = n≥1 1 n P (X (i) n = (k 1 , . . . , k i-1 , -1, k i+1 , . . . , k d )) ,
which gives (2.4) from (2.2). If (B i ) holds, then by definition, individuals of type i in F are all leaves and hence, X i,j ≡ 0, for all j = i and X i,i n = -n, for all n ≥ 0, see (2.1). In this case, the conclusion follows immediately.

Let us now prove properties 1-3 of µ. First note that for all i = j, m ij = 0 if and only if mij = 0. Then let r ≥ 1, assume that µ i admits moments of order r and that there is

j = i such that m ij = E(X i,j 1 ) > 0. The variable τ (i)
1 is a stopping time in the filtration generated by X (i) to which the increasing random walk X i,j is adapted. Then by applying Theorem 5.4 in [START_REF] Gut | Stopped random walks. Limit theorems and applications[END_REF], we obtain that E((X i,j 1 ) r ) < ∞ and E((τ

(i) 1 ) r ) < ∞. In particular τ (i) 1 < ∞, a.s. Now by definition, the random walk (X i,i n ) can be written as X i,i n = Y i,i n -n, where (Y i,i n ) is an increasing random walk. Since Y i,i (τ (i) 1 ) = τ (i) 1 -1 and E((τ (i) 1 ) r ) < ∞, we have E |Y i,i (τ (i)
1 )| r < ∞ and by applying Theorem 5.4 in [START_REF] Gut | Stopped random walks. Limit theorems and applications[END_REF] again, we obtain that

E |Y i,i 1 | r < ∞, and hence E |X i,i 1 | r < ∞.
So we have proved that ν admits moments of order r. Then it follows from the definition of τ (i) 1

and from Lemma 3.1 in [START_REF] Kesten | Two renewal theorems for general random walks tending to infinity[END_REF] that E((τ

(i) 1 ) r ) < ∞ implies lim n→∞ X i,i n = -∞,
and hence m ii < 1, from the law of large numbers.

Conversely, if m ij = 0 for all j = i, then mij = 0 for all j = i and µ i is the Dirac mass at 0, so it admits moments of order r. Now assume that ν i admits moments of order r and m ii < 1. Then it follows directly from Lemma 3.1 in [START_REF] Kesten | Two renewal theorems for general random walks tending to infinity[END_REF] that E((τ

(i) 1 ) r ) < ∞. Moreover from Theorem 5.2 in [10], E(X i,j (τ (i) 1 ) r ) < ∞,
for all j = i, which means that µ i admits moments of order r. If ν i admits moments of order 1 and m ii < 1, then it follows from the optional stopping theorem applied to the martingale (X i,j n -nE(X i,j 1 )), that E(X i,i (τ

(i) 1 )) = -1 = E(X i,i 1 )E(τ (i) 1 ) = (m ii -1)E(τ (i) 1
), and when i = j, E(X i,j (τ

(i) 1 )) = E(X i,j
1 )E(τ

(i) 1 ) = m ij
1-m ii and part 1 is proved. If M is irreducible, then for all i, there is j = i such that mij > 0. From part 1., ν i admits moments of order 1 and m ii < 1, for all i. In this case,

M + ∆ 2 = ∆ 1 M , where ∆ 1 = diag( 1 1 -m ii ) and ∆ 2 = diag( m ii 1 -m ii ) ,
and we derive from this identity that M is irreducible. Conversely if M is irreducible, then for all i, there is j = i such that m ij > 0 and hence mij > 0. Since by assumption, mij < ∞, for all i, j, then from part 1., m ii < 1, and M + ∆ 2 = ∆ 1 M holds. We derive from this identity that M is irreducible. Now if M is primitive, then it is irreducible and as before,

m ii < 1 for all i ∈ [d]. Moreover, M = (I -diag(m ii ))M + diag(m ii ) .
Therefore M is primitive. The converse cannot be true since there are nonnegative, irreducible matrices whose main diagonal is zero and which are not primitive. We can find distributions ν such that it is the case for M and hence for (I -diag(m ii ))M . If m ii > 0, for all i, then it follows from general theory of nonnegative matrices that M = (I -diag(m ii ))M + diag(m ii ) becomes primitive, see [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF]. Let us now prove 3. Recall that by definition, since M is primitive, µ i admits moments of order 1 for all i ∈ [d]. Then from the same arguments as in part 2.,

M = (I -diag(m ii ))M + diag(m ii ) and m ii < 1 for all i ∈ [d].
Assume that M is surpercritical, then there is a positive vector x such that M x > x. Therefore, (I -diag(m ii ))M x > (I -diag(m ii ))x and since m ii < 1, we obtain M x > x. Hence M is supercritical. Conversely, assume that M is supercritical. Then there is a positive vector x such that M x > x, so that

M x = (I -diag(m ii ))M x + diag(m ii )x > (I - diag(m ii ))x + diag(m ii )x = x and thus M is supercritical. Then the identity M = (I -diag(m ii ))M + diag(m ii ) allows us to derive that M is critical if and only if this is the case for M . Finally assume that m ii > 1 for some i ∈ [d]. If m ij = 0, for all j = i, then it is clear that individuals of type i in F are leaves. If m ij > 0, for some j ∈ [d]
, then since clusters of type i are supercritical, some of them have infinitely many children with positive probability. Conditionally to this event, such a cluster produces almost surely infinitely many children of type j, which is equivalent to say that individuals of type i in F give birth to an infinite number of children of type j with positive probability. 2

Let us now consider a multitype branching forest F with progeny distribution ν, with a finite number of trees and let

Z n = (Z (1) n , . . . , Z (d) n ), n ≥ 0 be the associated branching process, that is for each i ∈ [d], Z (i)
n is the total number of individuals of type i present in F at generation n. For x = (x 1 , . . . , x d ) ∈ Z d + , we denote by P x the law on (Ω, F) under which F is issued from x. In particular, P x (Z 0 = x) = 1. Then the next result gives the law of the total number of mutations in the forest F , that is the number of mutations up to the last generation whose rank is the extinction time, T := inf{n : Z n = 0}. For i, j ∈ [d], let us denote by M i the total number of mutations of type i in F , up to time T and by M ij the total number of mutations of type j produced by individuals of type i. In particular, M ii = 0 and M i and M ij satisfy the relations

M j = d i=1 M ij , j ∈ [d].
Note that if ν is primitive and supercritical, then P x (T = ∞) > 0 for all x ∈ Z d + , so that under P x , M i and M ij are infinite with positive probability, for some i, j ∈ [d]. We also emphasize that M i and M ij are not functionals of the branching process (Z n ).

Corollary 2.2. Assume that

(A i ) or (B i ) holds for all i ∈ [d]. Then for all integers x i , n i , k ij , i, j ∈ [d], such that x i ≥ 0, n i = -k ii , for i = j, k ij ≥ 0, and for all j ∈ [d], n j = x j + i =j k ij , P x (M 1 = n 1 -x 1 , . . . , M d = n d -x d , M ij = k ij , ∀i = j) = det(K) n1 . . . nd d i=1 µ * n i i (k i1 , . . . , k i(i-1) , 0, k i(i+1),...,k id ),
where µ i is defined in Theorem 2.1 and

µ * 0 i = δ 0 , ni = n i ∨ 1, K is the matrix (-k ij ) i,j
to which we removed the line i and the column i for all i such that n i = 0.

Proof.

This result is a direct consequence of Theorem 1.2 in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF] and Theorem 2.1 applied to the forest of mutations F . Indeed, it suffices to note that x i +M i corresponds to the total number of individuals of type i in F . Note however that Theorem 1.2 in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF] is proved only in the case where ν is primitive and (sub)critical. But using the coding which is presented in Section 2.1 and appyling Lemma 1 in [START_REF] Chaumont | Breadth first search coding of multitype forests with application to Lamperti representation[END_REF], we can check that it is still valid in the general case by following along the lines the proof which is given in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF]. 2

If for some i ∈ [d], none of the conditions (A i ) and (B i ) holds, then the definition of the vector of mutation sizes (M 1 , . . . , M d ) still makes sense. In this case, it is possible to obtain its law by extending Theorem 2.1 to branching forests whose progeny laws give mass to infinity. Note also that Corollary 2.2 can be considered as an extension of Theorem 1 in [START_REF] Bertoin | The structure of the allelic partition of the total population for Galton-Watson processes with neutral mutations[END_REF], where a similar formula is given in the neutral case.

We now turn our attention to the asymptotic behaviour of the number of mutations, when the total population is growing to infinity. Our first result is concerned with the critical case and is a direct consequence of Proposition 2 in [START_REF] Penisson | Various ways of conditioning multitype Galton-Watson processes[END_REF] and Theorem 2.1. If M is primitive, then we denote by u and v the unique right and left positive eigenvectors of M which are associated to the eigenvalue 1 and normalized by u.1 = u.v = 1. Recall that, for a multitype branching forest F , when no confusion is possible, N i denotes the total population of type i in F and M i denotes the total number of mutations of type i in F . Note also that when ν is primitive and critical, then (A i ) necessarily holds for all i ∈ [d], so that from Theorem 2.1, the forest of mutations F associated to F is a branching forest with progeny distribution µ defined by (2.4).

Corollary 2.3. Let F be a branching forest with a non singular, primitive and critical progeny distribution ν. Assume that for all i ∈ [d], µ i admits moments of order d + 1. If moreover M is primitive and the covariance matrices Σ i , Σ i of ν i and µ i , respectively are positive definite. Then m ii < 1, for all i ∈ [d] and there are constants

C 1 , C 2 > 0 such that for all x 0 ∈ Z d + , lim n→∞ n d/2+1 P x 0 (M i = n(1 -m ii )v i , i ∈ [d]) = C 1 x 0 .u , lim n→∞ n d+1 P x 0 (M i = n(1 -m ii )v i , N i = nv i , i ∈ [d]) = C 2 x 0 .u .
Proof. Since by assumption, M is primitive, then for all i, there is j = i such that mij > 0, and hence m ij > 0. Therefore, from part 1. of Theorem 2.1, m ii < 1, for all i. Moreover, from our assumptions and part 3. of Theorem 2.1, µ is critical. Besides, it is plain that M is non singular. Then conditions of Proposition 2 in [START_REF] Penisson | Various ways of conditioning multitype Galton-Watson processes[END_REF] are satisfied for the multitype branching process associated to F and the first assertion follows with ū and v, the normalized, positive right and left eigenvectors of M associated to the eigenvalue 1. Then recall from the proof of part 3. of Theorem 2.1 that M = (I -diag(m ii ))M + diag(m ii ). We derive from this identity that ū = u and v = cv(I -diag(m ii )), where c = u • v(I -diag(m ii )) -1 and the first assertion follows.

The proof of the second assertion follows the same lines as the proof of Proposition 2 in [START_REF] Penisson | Various ways of conditioning multitype Galton-Watson processes[END_REF]. In this case, since the number of mutations is taken into account together with the total number of individuals, a 2d-dimensional random walk is involved in the proof, which explains that the rate of convergence in now d + 1. 2 Note that the constants C 1 and C 2 can be made explicit in terms of the distributions ν and µ by properly exploiting the proof of Proposition 2 in [START_REF] Penisson | Various ways of conditioning multitype Galton-Watson processes[END_REF].

Through the next result we focus on the asymptotic behaviour of the number of mutations in a branching forest when the initial number of individuals x = (x 1 , . . . , x d ) tends to infinity along some given direction.

Theorem 2.4. Let F (x) be any family of multitype branching forests defined on the space (Ω, F, P ), indexed by x ∈ Z d + and such that for each x, F (x) has progeny distribution ν and is issued from x. For i ∈ [d], let N i (x) (resp. M i (x)) be the total number of individuals (resp. of mutations) of type i in F (x). Assume that ν is primitive and

let w ∈ Z d + \ {0}. 1. If ν is critical, then lim n→∞ N i (nw) n = ∞ and lim n→∞ M i (nw) N i (nw) = 1 -m ii , in probability. 2. If ν is subcritical, then lim n→∞ N i (nw) n = c i (w) and lim n→∞ M i (nw) n = w i + (1 -m ii )c i (w) , in probability, where c i (w) := d k=1 w k (I -M ) -1 ki . In any case, m ii < 1, for all i ∈ [d].
Proof. In order to prove our result, it suffices to construct some particular family of forests F (x), such that for each x, F (x) has progeny distribution ν and is issued from x ∈ Z d + , and to show that the limits in the statement hold.

Recall the coding of multitype branching forests which is presented at the end of Section 2.2 and let X (i) = {X i,j , j ∈ [d]} be d independent random walks whose respective step distributions are νi , i ∈ [d] defined in (2.2). Then for each x ∈ Z d + , we construct a forest F (x) such that F (x) is encoded by the random walks X (i) , i ∈ [d] and contains exactly x i trees whose root is of type i. This construction is possible in the primitive, (sub)critical case, thanks to part 3. of Theorem 3.1 in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF].

Then N i (x) and X (i) , i ∈ [d], satisfy identity (2.3). Moreover, for k = i, the number of mutations of type i issued from an individual of type k is X k,i (N k (x)), so that the total number of mutations of type i is

M i (x) = k =i X k,i (N k (x)) = -x i -X i,i (N i (x)) .
We derive from Lemma 2.2 in [START_REF] Chaumont | Coding multitype forests: application to the law of the total population of branching forests[END_REF], that if

x 1 , x 2 ∈ Z d + are such that x 1 ≤ x 2 , then the couple of random variables (N i (x 2 ) -N i (x 1 ), X i,i (N i (x 2 ) -X i,i (N i (x 1 ))) is inde- pendent of process ((N i (x), X i,i (N i (x))), x ≤ x 1 )
and has the same law as (N i (x 2x 1 ), X i,i (N i (x 2 -x 1 )). Therefore, for any w ∈ Z d + , ((N i (nw), X i,i (N i (nw)), n ≥ 0) is a bivariate random walk whose step distribution is the law of (N i (w), X i,i (N i (w)).

Let Z = (Z (1) , . . . , Z (d) ) be the branching process associated to F (w). Then by definition of N i (w), we have

N i (w) = ∞ n=0 Z (i) n . But E w (Z n ) = wM n , so that E w (Z (j) n ) = d i=1 w i m (n)
ij and since ν is primitive, we have from Frobenius Theorem for primitive matrices, m [START_REF] Athreya | Branching Processes[END_REF]. So we have proved that E(N i (w)) < ∞ if and only if ν is subcritical. Moreover, if ν is subcritical, then I -M is invertible and it follows from the above expressions that

(n) ij ∼ u i v j ρ n , see Theorem 1, Section V.2 in
E(N i (w)) = d i=1 w i (I -M ) -1 ij .
Then assertions 1. and 2. follow directly from the law of large numbers.

Finally, since ν and µ are primitive, by definition, they admit moments of order 1 and we derive from part 1. of Theorem 2.1 that m ii < 1, for all i ∈ Then this individual dies at the same time it gives birth. We emphasize that in this model, the probability for the population to become extinct does not depend on the rates λ i . This model is represented as a plane forest with edge lengths, see Figure 3. (In each sibling, we rank individuals of type 1 to the left, then individuals of type 2, and so on.) Such a forest will be called a multitype branching forest with edge lengths issued from x = (x 1 , . . . , x d ), with progeny distribution ν := (ν 1 , . . . , ν d ) and reproduction rates (λ 1 , . . . , λ d ). By construction, its discrete time skeleton is a multitype branching (plane) forest, as defined in the previous section, with progeny distribution ν, which is independent from the edge lengths. Edge lengths are independent between themselves and the length of an edge issued from a vertex of type i follows an exponential distribution with parameter λ i . We emphasize that the total number of individuals and the total number of mutations in a multitype branching forest with edge lengths are the same as in its discrete skeleton. Hence, the results of the previous section can be applied in the present setting.

Given a branching forest with edge lengths, as defined above, we denote by Z = (Z (1) , . . . , Z (d) ) the corresponding multitype branching process, that is for t ≥ 0 and i ∈

[d], Z (i)
t is the number of individuals of type i at time t in the population. (Since no confusion is possible, for the branching process we have kept the same notation as in discrete time.) The process Z is a Z d + -valued continuous time Markov process which satisfies the branching property, i.e., for λ ∈ R d + , t ≥ 0 and x, y ∈ Z d + ,

E x+y (e -λZt ) = E x (e -λZt )E y (e -λZt ) , where P x is the law under which the forest is issued from x. In particular, Z 0 = x, P xa.s. The process Z actually contains much less information than the original branching forest. In order to preserve the essential part of this information, we need to decompose Z as in the following definition. Definition 3.1. For i = j, we denote by Z i,j t the total number of individuals of type j whose parent has type i and who were born before time t. For i = j, the definition of Z i,i t is the same, except that to this number we add the number of individuals of type i at time 0 and we subtract the number of individuals of type i who died before time t.

The processes Z i,j whose definition should be clear from the example given in Figure 3 will play a crucial role in our continuous time model. A more formal definition can be found in Section 4.2 of [START_REF] Chaumont | Breadth first search coding of multitype forests with application to Lamperti representation[END_REF]. The interest of these processes is the following straightforward decomposition of the branching process Z = (Z (1) , . . . , Z (d) ):

(3.6) Z (j) t = d i=1 Z i,j t , j ∈ [d] .
Our model bears on a Lamperti type representation of these processes. According to Lamperti representation, any one dimensional branching process can be expressed as a Lévy process time changed by some integral functional. In this subsection, we will recall from [START_REF] Chaumont | Breadth first search coding of multitype forests with application to Lamperti representation[END_REF] the extension of this transformation to multitype, continuous time, discrete valued branching processes. The latter involves time changed multidimensional compound Poisson processes which we now introduce.

Since our models of evolution are only concerned with mutations, individuals of type i having exactly one child of type i do not present any interest. Hence we can assume without loss of generality that

ν i (e i ) = 0, for all i ∈ [d].
Then let X = (X (1) , . . . , X (d) ), where X (i) , i ∈ [d] are d independent Z d -valued compound Poisson processes. We assume that X (i) 0 = 0 and that X (i) has rate λ i and jump distribution νi which has been defined in (2.2). In particular, with the notation X (i) = (X i,1 , . . . , X i,d ), the process X i,i is a Z-valued, downward skip free, compound Poisson process, i.e. ∆X i,i t = X i,i t -X i,i t-≥ -1, t ≥ 0, with X 0-= 0 and for all i = j, the process X i,j is an increasing compound Poisson process. We emphasize that in this definition, some of the processes X i,j , i, j ∈ [d] can be identically equal to 0.

The following extension of the Lamperti representation to multitype branching processes can be found in [START_REF] Chaumont | Breadth first search coding of multitype forests with application to Lamperti representation[END_REF], see also [START_REF] Caballero | Affine processes on R n + × R n and multiparameter time changes[END_REF] for the case of continuous state multitype branching processes. Theorem 3.2. Let us consider a multitype branching forest with edge lengths issued from x = (x 1 , . . . , x d ) ∈ Z d + , with progeny distribution ν := (ν 1 , . . . , ν d ) and reproduction rates (λ 1 , . . . , λ d ). Then the processes Z i,j , i, j ∈ [d] introduced in Definition 3.1 admit the following representation:

(3.7) Z i,j t =      X i,j t 0 Z (i) s ds , t ≥ 0 , if i = j, x i + X i,i t 0 Z (i) s ds , t ≥ 0 , if i = j,
where the processes,

X (i) = (X i,1 , X i,2 , . . . , X i,d ) , i = 1, . . . , d ,
are independent Z d + valued compound Poisson processes, with jump distribution (ν 1 , . . . , νd ) and rates (λ 1 , . . . , λ d ). In particular from (3.6) and (3.7), the multitype branching process Z admits the following representation,

(3.8) (Z (1)
t , . . . , Z

t ) = x + d i=1 X i,1 t 0 Z (d)
, t ≥ 0 .

3.2.

Further results on asymptotics of mutations. For i ∈ [d] and t ≥ 0, we will denote by M i,t the total number of mutations of type i which occured up to time t. The definition of this quantity is illustrated on Figure 3. Let us also define a cluster of type i as the subtree corresponding to the descendence of type i of an individual of type i which is either a root or an individual whose parent as a type different from i. Then x i +M i,t corresponds to the number of clusters of type i in the forest truncated at time t.

In Proposition 3.4, we describe the asymptotic behaviour of M i,t , as t tends to ∞ in the case where the progeny distribution ν is primitive and supercritical. To this aim, we will need the joint representation of M i,t together with the number Z (i) t of individuals of type i at time t which is presented in Proposition 3.3.

Proposition 3.3. Recall from Section 3.1 the definition of the compound Poisson pro- cesses X i,j , i, j ∈ [d]. Then for any x = (x 1 , . . . , x d ) ∈ Z d + , under P x , the stochastic process Z (i) t , M i,t fulfills the following representation, Z (i) t , M i,t = x i + d k=1 X k,i t 0 Z (k) u du , d k=1,k =i X k,i t 0 Z (k) u du , t ≥ 0.
Proof. This result is a direct consequence of the representation which is recalled in Theorem 3.2. Indeed, recall from Section 3.1 the definition of Z i,j , then the number of mutations of type i up to time t is

M i,t = k =i Z k,i t .
The result follows from identity (3.

= 6, Z (1) t
= 3, Z 1,1 t = -2, Z 1,2 t = 5, Z 2,1 t = 8, Z 2,2 t = -2, and M 1,t = 8, M 2,t = 5. (2) t
Let's us now turn to the limiting behavior of M i,t , as t tends to infinity. The next result is concerned with the case where ν is primitive and supercritical. It allows us to evaluate the number of mutations which occured up to time t (or equivalently the number of clusters in the forest truncated at time t), when t is large.

Let us define the matrix A = Λ(M -I), where Λ = diag(λ i ). If M is primitive, then so is A and it follows from Perron-Frobenius theory that the eigenvalues ρ i , i ∈ [d] of A can be arranged so that ρ 1 > Re(ρ 2 ) ≥ • • • ≥ Re(ρ d ). Moreover, ν is subcritical, critical or supercritical according as ρ 1 < 0, ρ 1 = 0 or ρ 1 > 0. Then a well known result due to [START_REF] Athreya | Some results on multitype continuous time Markov branching processes[END_REF], see also Theorem 2, p. 206 in [START_REF] Athreya | Branching Processes[END_REF] asserts that when ν is non singular and primitive, there exists a nonnegative random variable W such that for all i ∈

[d], lim t→∞ e -ρ 1 t Z (i) t = v i W , a.s., (3.9)
where v i is the i-th coordinate of the normalized left eigenvector associated with ρ 1 . Proposition 3.4. Assume that ν is non singular, primitive and supercritical. Then

for all i ∈ [d], lim t→∞ e -ρ 1 t M i,t = K i W, a.s.,
where

K i = v i (1 + (1 -m ii )(λ i ρ 1 ) -1 ).
Proof. We derive from Proposition 3.3 that,

Z (i) t -M i,t = X i,i t 0 Z (i)
u du , a.s.

On the other hand, in the supercritical case, ρ 1 is strictly positive. Hence it follows from (3.9) that t 0

Z (i) u du ∼ ρ -1 1 W v i e ρ 1 t , a.s., as t → ∞.
Then the desired result is a consequence of the latter equivalence and the law of large numbers applied to the compound Poisson process X i,i . 2

Under conditions of Propositon 3.4, assume moreover that for some i ∈

[d], K i is positive, that is m ii < 1 + λ i ρ 1 ,
and that for some j, P e j (W > 0) = 1. Then using Proposition 3.4, we can compare the asymptotic behaviour of the number of mutations prior to t with this of Z (i) t , under P e j , that is

(3.10) M i,t ∼ K i Z (i)
t , P e j -a.s., as t → ∞. Regarding the condition P e j (W > 0) = 1, note that Theorem 2, p. 206 in [START_REF] Athreya | Branching Processes[END_REF] also asserts that P e k (W > 0) > 0, for some (hence for all) k ∈ [d], if and only if

(3.11) E(ξ ij log ξ ij ) < ∞, for all i, j ∈ [d],
where (ξ i1 , . . . , ξ id ) is a random vector with law ν i . Moreover, 1 -P e k (W > 0) corresponds to the probability of extinction, when the forest is issued from e k .

3.3. Emergence times of mutations. In this section, we shall assume that mutations are not reversible, that is for all i = 1, . . . , d -1, individuals of type i can only have children of type i or i+1. In particular ν is not irreducible. Moreover when giving birth, individuals of type i = 1, . . . , d -1 have at least one child of type i with probability one, and have children of type i + 1 with positive probability. These conditions can be explicited in terms of the progeny distribution ν i as follows (3.12)

   ν i (k) > 0 ⇒ k j = 0, for j / ∈ {i, i + 1}, k∈Z d + :k i =0 ν i (k) = 0 and k∈Z d + :k i+1 =0 ν i (k) < 1.
We are interested in the waiting time until an individual of type i first emerges in the population, that is

τ i := inf{t ≥ 0 : Z (i) t ≥ 1} .
The problem of determining a general expression for the law of τ i is quite challenging. As far as we know, there is no explicit expression for this law in terms of the progeny distribution and the reproduction rates. Various results in this direction can be found in [START_REF] Serra | On the Waiting Time to Escape[END_REF], [START_REF] Serra | Dynamics of escape mutants[END_REF], [START_REF] Durrett | Evolution of resistance and progression to disease during clonal expansion of cancer[END_REF] and [START_REF] Alexander | Conditional distributions and waiting times in Multitype branching processes[END_REF] for instance. Most of them provide approximations of this law, using martingale convergence theorems [START_REF] Durrett | Evolution of resistance and progression to disease during clonal expansion of cancer[END_REF] or through numerical methods for the inversion of the generating function [START_REF] Alexander | Conditional distributions and waiting times in Multitype branching processes[END_REF]. In Proposition 3.5 we first give a relationship between the successive emergence times τ 2 , τ 3 , . . . in terms of the underlying compound Poisson process in the Lamperti representation of Z. We also characterize the joint law under P e i-1 of the time τ i and the number of individuals of type i -1 at this time. In Theorem 3.7 we derive an approximation of the time τ i , under P e 1 , as the mutation rate of type k increases faster than that of type k -1, for all k = 3, . . . , i. Then in Corollary 3.8 we focus on a case where these law can be explicited.

In the following developments, we use the notation of Section 3.1 from which we recall the Lamperti representation of the multitype branching process Z = (Z (1) , . . . , Z (d) ) in terms of the compound Poisson processes X (i) . Let us also introduce a few more notation. For i, j ∈ [d], we denote by λ i,j the parameter of the compound Poisson process X i,j , that is λ i,j := λ i 1 -k∈Z d + :k j =0 νi (k) . Note that from our assumptions (3.12), for all i = 1, . . . , d -1, λ i,i+1 > 0 and for j / ∈ {i, i + 1}, λ i,j = 0, that is X i,j is identically equal to 0. In particular, λ i = λ i,i + λ i,i+1 , for i ≤ d -1 and λ d = λ d,d . The parameter λ i,i+1 will be call the mutation rate of type i + 1. For i ≥ 2, let γ i := inf{t : X i-1,i t ≥ 1} be the time of the first jump by the process X i-1,i and note that this time is exponentially distributed with parameter λ i-1,i . Proposition 3.5. Assume that (3.12) holds and define Z 0,1 as the process identically equal to 1 and set τ 1 = 0.

1. For i = 2, . . . , d, the emergence time τ i of type i admits the following representation under P e 1 , (3.13)

τ i = τ i-1 + γ i 0 1 X i-1,i-1 s + Z i-2,i-1 κ i-1 (s)
ds, P e 1 -a.s., where κ i-1 is the right continuous inverse of the functional t → t 0 Z (i-1) s ds, i.e. κ i-1 (t) = inf{s > 0 :

s 0 Z (i-1) u du > t}. 2.
Under P e i-1 , the joint law of the emergence time τ i of type i together with the number of individuals of type i -1 in the population at time τ i admits the following representation, θ k > t , for all t > 0 .

(3.14) (τ i , Z (i-1) τ i ) (d) = γ i 0 ds 1 + X i-1,i-1 s , 1 + X i-1,i-1 γ i . 3. Let us define θ k = γ k 0 1 X k-
Proof. Since X i,j is identically equal to 0 whenever j / ∈ {i, i + 1}, then under P e 1 , the representation (3.8) admits the simpler form (3.16)

(Z (1) t , . . . , Z (d) t ) = e 1 + X 1,1 t 0 Z (1) s ds , X 2,2 t 0 Z (2) s ds + X 1,2 t 0 Z (1) s ds , . . . , X d,d t 0 Z (d) s ds + X d-1,d t 0 Z (d-1) s ds .
In particular, for i = 2, . . . , d,

Z (i) t = X i,i t 0 Z (i) s ds + X i-1,i t 0 Z (i-1) s ds , t ≥ 0.
Since X i,i 0 = 0, for i ≥ 2, we see that the time τ i corresponds to the first hitting time of level 1 by the process t →

X i-1,i t 0 Z (i-1) s ds , that is (3.17) τ i = κ i-1 (γ i ),
where γ i has been defined as the time of the first jump of the process X i-1,i . For t such that κ i-1 (t) < ∞, we have t =

κ i-1 (t) 0 Z (i-1) s ds, so that dt = Z (i-1) κ i-1 (t) dκ i-1 (t), and
since κ i-1 (0) = τ i-1 , we obtain κ i-1 (t) = τ i-1 + t 0 ds Z (i-1) κ i-1 (s) = τ i-1 + t 0 ds X i-1,i-1 s + X i-2,i-1 κ i-1 (s) 0 Z (i-2) u du . (3.18)
The latter identity together with (3.17) prove identity (3.13).

The second part of the proposition is easily derived from the same arguments. More specifically, it follows from (3.17) and the following identities

Z (i-1) t = 1 + X i-1,i-1 t 0 Z (i-1) s ds and κ i-1 (t) = t 0 ds 1 + X i-1,i-1 s , t ≥ 0 , which hold P e i-1 -a.s.
Independence between the variables θ k , k ≥ 2 is a direct consequence of the independence between the processes X (i) , i ∈ [d]. We derive from the representation of τ i in part 1. of this proposition that (3.19) Note that the law of θ k or equivalently, the law of τ k under P e k-1 can be made explicit in some instances through its Laplace transform, see Corollary 3.8 below.

τ i = i k=2 γ k 0 1 X k-1,k-1 s + X k-2,k-1 κ k-1 (s) 0 Z (k-2) u du ds, a.s. Note that since κ k-1 (0) = τ k-1 , then from (3.17), for all k ≥ 2, κ k-1 (0) 0 Z (k-2) u du = γ k-1 , so that by definition of γ k-1 , (3.20) X k-2,k-1 κ k-1 (0) 0 Z (k-2) u du = X k-2,k-1 γ k-1 ≥ 1, a.s. Besides, since s → X k-2,k-1 κ k-1 (
For the remainder of this section we will assume moreover that at each mutation, individuals of type i do not give birth to more than one child of type i + 1 in a same litter. More specifically, assumptions (3.12) are replaced by,

(3.21)    ν i (k) > 0 ⇒ k i+1 = 0 or 1 and k j = 0, for j / ∈ {i, i + 1}, k∈Z d + :k i =0 ν i (k) = 0 and k∈Z d + :k i+1 =0 ν i (k) < 1 .
In particular, under these assumptions, the process X i,i+1 is a standard Poisson process. Then we will need the next lemma in order to derive our main result on the estimation of the time τ i , as the mutation rates λ k-1,k , k = 2, . . . , d grow faster. Lemma 3.6. Assume that (3.21) holds, let k ≥ 3 and fix λ 1,2 > 0, then

P e 1 X k-2,k-1 κ k-1 (γ k ) 0 Z (k-2) u du = 1 -→ 1 , as λ n-2,n-1 /λ n-1,n → 0, for n = 3, . . . , k. (3.22)
increase very fast. It would interesting to study the asymptotic behavior of τ i , when λ k,k λ k+1,k+1 → 0, that is when the intrinsic reproduction rates increase very fast. This assumption also fits to the model of cancer, since mutations are always more sensitive to proliferate.

λ k-2,k-1 λ k-1,k → 0, for k = 3, . . . , i.
Besides, the expectation of τ i fulfills the following approximation:

E e 1 (τ i ) ∼ i k=2 E(θ k ) , as λ k-2,k-1 λ k-1,k → 0, for k = 3, . . . , i. Proof. Since s → X k-2,k-1 κ k-1 (s) 0 Z (k-2) u du
are increasing processes, then from (3.20),

P e 1 -almost surely on the set {X k-2,k-1 κ k-1 (γ k ) 0 Z (k-2) u du = 1}, we have γ k 0 1 X k-1,k-1 s + X k-2,k-1 κ k-1 (s) 0 Z (k-2) u du ds = γ k 0 1 X k-1,k-1 s + 1 ds .
Hence it follows from Lemma 3.6 that for fixed λ 1,2 > 0, as λ n-2,n-1 /λ n-1,n → 0, for all n = 3, . . . , k,

γ k 0 1 X k-1,k-1 s + 1 ds -1 γ k 0 1 X k-1,k-1 s + X k-2,k-1 κ k-1 (s) 0 Z (k-2) u du ds P -→ 1 ,
and the first part of the theorem is easily derived from this convergence and (3.13) (or equivalently (3.19)).

In order to prove the second part, let us first set But for any p, q ≥ 1, such that p -1 + q -1 = 1, we have from Holder inequality E(θ k 1 A c k ) ≤ E(θ p k ) 1/p P (A c k ) 1/q . Moreover, we clearly have E(θ p k ) 1/p ∼ 1/λ k-1,k , as λ k-1,k → ∞. Hence, (3.27) is satisfied thanks to Lemma 3.6. 2

H k := γ k 0 1 X k-1,k-1 s + X k-2,k-1
We end this section with an example where the distribution of τ i can be estimated a bit more specifically. We consider the case of binary fission with mutations, where each individual of type i can give birth to either two individuals of type i or one individual of type i and one individual of type i + 1. In particular, all jumps of Z i,i have size 1 and X i,i is a Poisson process with parameter λ i,i . Corollary 3.8. With the above assumtions, the law of τ i can be specified as follows.

1.

Under P e i-1 , the Laplace transform of τ i is expressed as,

E e i-1 (e -ατ i ) = λ i-1,i n≥0 λ n i-1,i-1 n k=0 (α k + • • • + α n + ᾱn+1 )
, α ≥ 0 , where α 0 = 0, α k = α k(k+1) and ᾱk = λ i-1 + α k , for k ≥ 1. 2. The expectation of τ i is given by E e i-1 (τ i ) = 1 λ i-1,i λ i-1,i-1 ln λ i-1 λ i-1,i . In particular, for fixed λ 1,2 > 0, under P e 1 , the expectation of τ i fulfills the following approximation:

E e 1 (τ i ) ∼ i k=2 λ -2
k-1,k , as Under P e i-1 , X i-1,i-1 is a standard Poisson process with parameter λ i-1,i-1 starting at 0. So if we denote by (J n ) n≥1 the sequence of jump times of X i-1,i-1 and set J 0 = 0, then developing the expression E e Then coming back to expression (3.28), we obtain with the convention that 0 k=1 = 0,

E e i-1 (e -ατ i ) = λ i-1,i n≥0 λ n i-1,i-1 0≤x 1 ≤•••≤x n+1 e -( ᾱn+1 x n+1 + n k=1 α k x k ) dx 1 . . . dx n+1 ,
where α 1 , . . . , α n , ᾱn+1 are defined in the satement. (Here we used the fact that λ i-1 = λ i-1,i + λ i-1,i-1 .) The computation of the integral is easily done.

Then using again part 2. of Proposition 3.5, we obtain the expectation of τ i under P e i-1 , after easy computations, E e i-1 (τ i ) = +∞ 0 dxλ i-1,i e -λ i-1,i x x 0 e -λ i-1,i-1 s k≥0 (λ i-1,i-1 s) k (k + 1)! ds = 1 λ i-1,i λ i-1,i-1 ln λ i-1 λ i-1,i .

  k , . . . the sequence of clusters of type i in (f , c f ). The forest f (i) := {t

[d]. 2 3.

 2 When continuous time is involved 3.1. The Lamperti representation. Let us now consider a d type population which is composed at time t = 0, of x i individuals of type i ∈ [d] and whose dynamics in continuous time behave according to a branching model. More specifically, at any time, all individuals in the population live, give birth and die independently of each other. Once it is born, any individual of type i ∈ [d] gives birth after an exponential time with parameter λ i > 0 to n j individuals of type j ∈ [d] with probability ν i (n 1 , . . . , n d ).

  , then inequality (3.15) is a direct consequence of identities (3.19) and (3.20). 2

Theorem 3 . 7 .

 37 Assume that (3.21) holds. Recall the definition of θ k in Proposition 3.5 and let us fix λ 1,2 > 0, then under P e 1 ,

1 κ

 1 ds and A k := {X k-2,k-Then from(3.13), E e 1 (τ i ) = i k=2 E e 1 (H k ), so it suffices to prove that for all k = 2, . . . , i,(3.26) E e 1 (H k ) ∼ E(θ k ), as λ n-2,n-1 /λ n-1,n → 0, for n = 3, . . . , k. Observe that E e 1 (H k ) = E(θ k 1 A k ) + E e 1 (H k 1 A c k ). Moreover, E e 1 (H k 1 A c k ) ≤ E(θ k 1 A c k ). Then to obtain (3.26), it is enough to prove that (3.27) E(θ k 1 A c k ) E(θ k )→ 0, as λ n-2,n-1 /λ n-1,n → 0, for n = 3, . . . , k.

λ k- 2 ,k- 1 λ

 21 k-1,k → 0, for k = 3, . . . , i.Proof. From part 2. of Proposition 3.5, for all β ≥ 0,E e i-1 (e -ατ i ) = E e i-1,i x dx . (3.28)

=

  e -(α+λ i-1,i-1 )x + n≥1 e -λ i-1,i-1 x (λ i-1,i-1 x) ) ) n! x n dx 1 . . . dx n = e -(α+λ i-1,i-1 )x + n≥1 λ n i-1,i-1 e -(λ i-1,i-1 + α n+1 ) dx 1 . . . dx n .

  2.1. Let F be any multitype branching forest with progeny distribution ν = (ν 1 , . . . , ν d ) and denote by F the associated forest of mutations. Assume that for all i ∈ [d], one of the conditions (A i ) or (B i ) holds. Then F is a multitype branching forest with progeny distribution µ = (µ 1 , . . . , µ d ) on S i := {k ∈ Z d + :

  7) in Theorem 3.2. 2

	t
	t=0
	Figure 3. A two type forest with edge lengths issued from x = (2, 2).
	Vertices of type 1 (resp. 2) are represented in black (resp. grey). At time
	t,
	Z

  +1 ds, for k ≥ 2.Then the random variables θ k , k ≥ 2 are independent and for i = 2, . . . , d,

		s	1,k-1
			i
	(3.15)	P e 1 (τ i > t) ≤ P
			k=2

We conclude from Theorem 3.7. 2

Proof. First set γ

(1)

= 2, P e 1 -a.s.) So from (3.17), we have showed that

k-1 } .

The event {τ k < τ

k-1 } means that before the first time when an individual of type k appears in the population, there has been only one birth of type k -1. From the Markov property applied at time τ k-1 , we have

The support in the integral of (3.24) is included in the set {z : z k-1 = 1}, so from (3.23), (3.24) and the Lebesgue theorem of dominated convergence, all we need to prove is

k-1 ) = 1, since in the first case Z k-2,k-1 is identically equal to 0, so that τ

k-1 = ∞, P z -a.s. and in the second case, τ k = 0, P z -a.s.) Let z be such that z k-1 = 1. Without loss of generality we can assume that z i ≥ 1, for i = 1, . . . k -2. For i = 1, . . . k -1, let us denote by U i the first time that the lineage of one of the z k-i initial individuals of type k -i gives birth to an individual of type k -i + 1. Then from the branching property, under P z , the r.v.'s U i are independent and from part 2. of Proposition 3.5, U i has the same law as

k-1 } , which imply the inequality,

k-1 ) . But when λ n-2,n-1 /λ n-1,n → 0, for n = 3, . . . , k, the parameter λ 1,2 > 0 being fixed, we necessarily have lim λ n-1,n = ∞, for n = 3, . . . , k. Hence γ k /γ k-1 converges in probability toward 0, 1/Y

(2)

), for n = 3, . . . , k -1 converge in probability toward +∞. Therefore, the left hand side of the above inequality tends to 1, which proves (3.25) and the lemma is proved. 2

In the following theorem, the assumption λ k-1,k λ k,k+1 → 0 is quite adapted to several biological models such as cancer growth, for instance. Indeed, cancer is often the result of a series of successive mutations, [START_REF] Iwasa | Population genetics of tumor suppressor genes[END_REF], [START_REF] Durrett | Evolution of resistance and progression to disease during clonal expansion of cancer[END_REF], [START_REF] Durrett | Population genetics of neutral mutations in exponentially growing cancer cell populations[END_REF]. Each new mutation is itself more unstable than the previous ones, and in particular, the successive mutation rates can


In [dMP13], L. de Moura and G. O. Passmore presented the Strategy Challenge to the SMT-community: "To build theoretical and practical tools allowing users to exert strategic control over core heuristic aspects of high-performance SMT solvers." They advocated for the adaptation of ideas of strategy prevalent in the Argonne and LCF theorem proving paradigms.

The same year, at the Tableaux conference, a system description was published of our Proof-Search factorY for Collaborative HEuristics, Psyche, entitled "A proof-search engine based on sequent calculus with an LCF-style architecture" [GL13]. Our highly modular system, designed for either interactive or automated theorem proving, provides a platform where users can experiment various techniques and strategies, to be implemented as plugins via an API with Psyche's kernel. Moreover, Psyche's framework features the possibility to call theory-specific decision procedures, so that the system can be used for SMTsolving: Indeed we showed in [FGLM13] that (an elementary version of) the DPLL(T ) scheme is bi-similar (i.e. isomorphic) to a simple proof-search strategy scheme, an instance of which we implemented as one of Psyche's plugins.

In [dMP13] is emphasized the key role that is played in SMT-solving by what is not specified by the DPLL(T ) scheme, such as heuristics. The ability to easily tune these and give users control over them, is a desirable feature to increase a solver's versatility. The Strategy Challenge suggests that this ambition deserves more than behaviour-controlling parameters to be set by the user; it may also require a more modular organisation of a solver's code, for instance to enhance the programmability of its heuristics. This is what we did with Psyche's API.

Of course, empowering users should not jeopardise the correctness of a solver's answers. But a modular architecture can in fact help, rather than jeopardise, the trust in a solver. And this is an extra motivation for the Strategy Challenge which adds on to [dMP13]. If the code that defines and controls the heuristics is separated from the code that implements only the high-level, non-deterministic, rule-based procedures such as DPLL(T ), then it is possible, using the right software abstractions, that only the latter kind of code needs to be trusted (or, even better, verified) in order to trust the whole solver. Regardless of the user's own fiddling of heuristics. This is where the LCF-inspired ideas come in: Psyche gives users full control in that they can write their own source code in the form of a plugin, as long as it implements a predefined interface. That interface requires their main function to output a value in the type of answers: type answer = Provable of statement*proof | NotProvable of statement

The only limitation, which guarantees correctness, is that this type of answers is private to Psyche's kernel: in order to inhabit it, the plugin cannot use its constructors and has no choice but to use the kernel's API, hoping that, from this interaction, a value in this type arises that can be given as the final output.

The approach we took in Psyche is to model this interaction as a slot machine, as each call to the kernel finishes with a value of the following recursive type which the plugin can exploit:

type output = Jackpot of answer | InsertCoin of info*(coin -> output) If the plugin is lucky, the kernel provides the answer the plugin is waiting for, otherwise it provides, besides some useful information, a function of type coin -> output which the plugin can use to 'put a new coin in the slot machine' instructing the kernel which rule to apply next (say in the DPLL(T ) scheme). And for this the plugin can summon all imaginable heuristics without jeopardising the correctness of an eventual answer.

This 'slot machine' form of kernel/plugin interaction, using a private type of answers to guarantee correctness, can be set up for any inference system in the sense of of [GRS04,Sha05], where the heuristics to be programmed as plugins are called reduction operations. Inference systems are used in particular as the basis of a unifying approach to describe theory combination techniques, so a slot machine implementation of them would also allow users to control how theories are combined in SMT-solving. This is one of Psyche's next developments.

The slot machine architecture used in Psyche gives users a control over strategies that is powerful, in that they have access to the full expressivity of the programming language to implement plugins, and fine-grained, in that they can decide up to the details of e.g. how and when unit propagation is performed. This "small step" strategic control over reasoning engines was left as future work in [dMP13]. Of course, the expressivity of a full programming language and the control of the small steps might be "too much" for the user who may not want (or know how) to harness so much power. In that case, one can provid intermediate modules between the kernel and what the user actually manipulates, for instance by nesting slot machines within slot machines: the deepest one interacts directly with the kernel on all heuristical aspects, while the top-level one interacts with the end-user on the much coarser heuristical aspects that [dMP13] mostly discusses. If need be, one may also offer a limited tactic-programming language for the end-user to specify how to insert coins into the slot machine, without exposing him to the programming language in which the prover is written.

Also, a key aspect of [dMP13] is the preservation of efficiency, especially for the small-step strategic control which may rely on the use of clever datastructures. More work needs to be done with Psyche to show that the slot machine architecture does not intrinsically yield unreasonable performance drops.

  Slot Machines: an approach to the Strategy Challenge in SMT solving (presentation only) Stéphane Graham-Lengrand SRI International -CNRS -École Polytechnique Abstract In this short introduction we briefly describe the relevance of Psyche's system description [GL13], attached, to the Strategy Challenge set by L. de Moura and G. O. Passmore for SMT solving.