Jonathan Passerat-Palmbach

Claude Mazel

Antoine Mahul

David R C Hill

Reliable Initialization of GPU-enabled Parallel Stochastic Simulations Using Mersenne Twister for Graphics Processors

Keywords: Parallel Stochastic Simulation, Random Numbers Generation, GP-GPU, CUDA, MTGP, TestU01

Parallel stochastic simulations tend to exploit more and more computing power and they are now also developed for General Purpose Graphics Process Units (GP-GPUs). Consequently, they need reliable random sources to feed their applications. We propose a survey of the current Pseudo Random Numbers Generators (PRNG) available on GPU. We give a particular focus to the recent Mersenne Twister for Graphics Processors (MTGP) that has just been released. Our work provides empirically checked statuses designed to initialize a particular configuration of this generator, in order to prevent any potential bias introduced by the parallelization of the PRNG.

INTRODUCTION

Over the past five years, simulationists have leant towards Graphics Process Units (GPUs) to compute the heavy tasks bound to their research activity. Stochastic simulations, especially Monte Carlo [START_REF] Gentle | Random Number Generation and Monte Carlo Methods (Statistics and Computing)[END_REF], are widely spread through many scientific communities. General-Purpose GPUs (GP-GPUs) with many cores are very interesting for this simulation technique. Depending on the domain in which this kind of simulation is used, results may suffer from a weak parallelization technique with some stochastic streams of poor quality. We have shown that for sensitive research domains, like nuclear medicine, parallel simulations cannot cope with biased results [START_REF] El Bitar | Fully 3d monte carlo image reconstruction in spect using functional regions[END_REF][START_REF] Reuillon | Rigorous distribution of stochastic simulations using the distme toolkit[END_REF].

Finding a fast and reliable Pseudo Random Number Generator (PRNG) to feed a sequential stochastic simulation is not a problem for many application domains since more than a decade. This issue has been tackled in many reference studies for CPU based PRNG [START_REF] L'ecuyer | Random numbers for simulation[END_REF]. [START_REF] Park | Random number generators: Good ones are hard to find[END_REF] raised the fact that one should consider the couple (application, PRNG) instead of limiting the study to the intrinsic qualities of the PRNG. For instance, a very good generator like Mersenne Twister (MT) (Matsumoto and Nishimura, 1998) should not be trusted for a cryptographic application and the initialization of such generators has to be done carefully. Indeed, for some years, this very nice generator was sensible to its initialization status. Even if we never have a universal generator, the MT family of generators [START_REF] Saito | Simd-oriented fast mersenne twister: a 128-bit pseudorandom number generator[END_REF], the WELL generators [START_REF] Panneton | Improved long-period generators based on linear reccurences modulo 2[END_REF] and some advanced Multiple Recursive pseudo-random numbers Generators (MRGs) from L'Ecuyer [START_REF] L'ecuyer | An object-oriented random-number package with many long streams and substreams[END_REF] give very good results when considered for parallel computing in a wide range of applications.

If some criteria have been gathered by Coddington [START_REF] Coddington | Random number generator for parallel computers[END_REF] for sequential and parallel simulation to characterize good PRNGs, it is often safer to refer to empirical testing software. Several libraries have been designed in this way, the oldest and most renowned include the statistical tests proposed by Knuth (Knuth, 1969), the DieHard testing suite designed by Marsaglia [START_REF] Marsaglia | The marsaglia random number cdrom, with the diehard battery of tests of randomness[END_REF] and the Brown's DieHarder suite [START_REF] Brown | Dieharder: A random number test suite[END_REF]. The TestU01 library provided by L'Ecuyer [START_REF] L'ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF] is the ultimate test software at the time of this paper writing. The latter contains a battery of tests called BigCrush, the most stringent set of tests that a PRNG can be faced nowadays.

The main problem we are still facing today is to ensure the correct behavior of the PRNG when distributed across a tight or large coupled computing architecture. The literature currently provides quite a few references about stochastic streams distribution on classical hardware architectures [START_REF] Mascagni | Some methods of parallel pseudorandom number generation[END_REF], [START_REF] Traore | The use of random number generation for stochastic distributed simulation: application to ecological modeling[END_REF], [START_REF] Bauke | Random numbers for large scale distributed monte carlo simulations[END_REF], [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF]. The set of references is even poorer when considering a GP-GPU environment. Indeed, restricted parallel hardware architectures like the Single Instruction Multiple Data (SIMD) family, which GPUs belong to, do highly impact the implementation of generators. In addition, we still have to select the best way to allocate random substreams to these manycore architectures.

In this paper, we bring our contribution to these problems:

• We survey the current PRNGs implementations and solutions available to distribute stochastic streams through a GPU architecture.

• The parameterized PRNGs family is proposed taking the example of MTGP [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF], the generator we find the most interesting in this context.

• We present the latter and discuss its quality with empirical benchmarks and statistical analysis. Reliable initialization data structures for a particular configuration of MTGP are also presented.

In this manner, we first intend to analyze the features of a particular generator designed for GPU hardware architectures: MTGP. The second purpose of this study is to give reliable parameters to initialize this PRNG, without introducing any potential bias in the parallel stochastic simulations based upon it.

PSEUDO RANDOM NUMBERS FOR GP-GPUs

Until recently, designing a PRNG for GPU-enabled platforms could be very tricky as it forced programmers to deal with graphics Application Programming Interfaces (APIs). Some implementations are presented in [START_REF] Sussman | Pseudorandom number generation on the gpu[END_REF]. The authors especially list the limitations of these GPU dedicated PRNGs due to the past weaknesses of the hardware. Limited output per thread or untruthful operations were part of the restrictions that made these PRNGs feeble for High Performance Computing (HPC) applications. Consequently, a common way to deal with random numbers on GPU was to generate them on CPU before transferring them on the graphics processor. This solution has to face the well-known bottleneck of data transfer between the CPU host and the GPU device. Even with nowadays PCI Express 16X running at 8GB/s, this remains a challenge for high performance applications.

Since 2008 and the recent advances from NVIDIA, new GPU software and hardware architectures offer the precision and speed needed for many HPC applications. Now, PRNGs can be directly implemented into the GPU. Recent works propose this new kind of generators. Langdon presents a minimal implementation of the standard Park Miller PRNG [START_REF] Park | Random number generators: Good ones are hard to find[END_REF] on a NVIDIA 8800 GTX GPU in its paper from 2008 [START_REF] Langdon | A fast high quality pseudo random number generator for graphics processing units[END_REF]. He announces a speed up of more than 40 compared to his Intel 2.40 GHz CPU. One year later, [START_REF] Langdon | A fast high quality pseudo random number generator for nvidia cuda[END_REF] increased again the speed of his application by four by using the new NVIDIA technology CUDA (Compute Unified Device Architecture) (NVIDIA, 2010) with a Tesla T10 GPU. Nevertheless, we do not advise the use of this old generator which has many known flows, though it was still in use until recently in some well distributed networking simulation software [START_REF] Entacher | Pitfalls when using parallel streams in omnet++ simulations[END_REF].

CUDA has been designed to allow developers to easily harness the computation power of GPUs. In his first implementation, Langdon had to deal with a complex and unadapted graphics API. With CUDA, developers can program GPUs without wasting their time making algorithms and their data fit into graphics dedicated data structures, such as pixels shaders. Furthermore, CUDA does not propose a new programming language but only some C extensions, making it easier to learn for C familiars. The CUDA appellation also refers to the name of the new NVIDIA GP-GPU architecture. This generation of graphic boards tries to fulfill the requirements noted in the conclusion of the previously cited [START_REF] Sussman | Pseudorandom number generation on the gpu[END_REF], with for instance an implementation of the IEEE 754 floating point numbers standard. The new generation of boards based upon the Fermi architecture is now proposing configurable L1 cache, ECC memory and a considerable increase of performance in double precision, while owning twice as much cores as the antecedently mentioned T10 processor.

Although these highly parallel devices bring much more peak performances than CPUs, they must be carefully programmed to deliver the expected power. In fact, GPU architectures combine a manycore approach with SIMD vector cores. As vector processors do, GPU-enabled algorithms need to repeat the same operation on different data to correctly exploit the device. This is the main reason of the recent dedicated PRNGs proposals. In 2006, [START_REF] Saito | Simd-oriented fast mersenne twister: a 128-bit pseudorandom number generator[END_REF] proposed an SIMD version of the famous 'Mersenne Twister' called SIMD-oriented Fast Mersenne Twister (SFMT). Although this algorithm can be used on regular CPUs or on SIMD enabled CPUs (using either SSE or AltiVec vector instructions), it cannot be directly transposed to a GPU architecture. Most PRNGs have to be rethought from scratch to leverage GPUs characteristics, once again, we always have to take into account the target application. In the case of a CUDA implementation, these couples are well surveyed in [START_REF] Howes | GPU Gems 3, chapter 37 -Efficient Random Number Generation and Application Using CUDA[END_REF].

Given that CUDA defines software levels that map the device architecture, PRNGs implemented using this technology can be organized at one of the following scopes, corresponding to the main elements of the CUDA framework: a thread, a block of threads or a kernel (the program running on the whole GPU). All these approaches have been studied in the literature. In [START_REF] Zhmurov | Efficient pseudo-random number generators for biomolecular simulations on graphics processors[END_REF], authors present three basic generation algorithms working either with a single instance of the PRNG for the kernel or with an instance per thread. The three algorithms exposed are quite basic: Ran2, Hybrid Taus and a Lagged Fibonacci generator. In the same way, [START_REF] Langdon | A fast high quality pseudo random number generator for nvidia cuda[END_REF] chooses to generate a number per thread in its GPU version of the Park-Miller algorithm. The last strategy is proposed in [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF] where a new variant of the MT algorithm spreads independent PRNGs through each thread block, thanks to an algorithm known as Dynamic Creator (DC) [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF] that we will detail later.

Beyond the nature of the PRNG algorithm, we prefer to focus on the scope chosen for each implementation. Indeed, we formerly insisted on the need to consider the target application and the PRNG as a pair. Obviously, new PRNG algorithms have to take advantage of GPU intrinsic properties such as heterogeneous memories, or thread organization. The former highly impacts the PRNG performance. Considering the approach using a generator per thread approach, an internal state array has to be saved in each thread. CUDA related works like [START_REF] Kirk | Programming Massively Parallel Processors[END_REF] specify that arrays declared for a thread are stored in the local memory, implemented in RAM. Equivalently, with a PRNG for all thread approach, the global memory is solicited to store the state of the PRNG. Each thread draws a number and updates its component of the state in global memory, implemented in RAM too. These two approaches make a heavy use of global memory, which has the advantage to be persistent across kernel launches within the same application. Yet, this RAM area is quite slow, it implies a 400 to 800 clock cycles latency because it is not cached (NVIDIA, 2010). So, even if the global memory storage is compulsory to save the PRNG state between two kernel calls, one can use the shared memory, reachable by every thread within a block, to manipulate PRNG data. Indeed, it is implemented on-chip and is consequently as fast as registers. A good example of this choice is the "to be published" paper introducing MTGP [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF].

In our opinion, a good GPU PRNG should employ shared memory. A PRNG per block approach seems to be the most appropriate way to implement a source of randomness, first, because it exploits the quickest memory, second, for the sake of applications upgradability. Since hardware architectures evolve very quickly, we cannot afford to rethink algorithms every time the memory amount or number of threads available doubles. So, fixing a block of threads grained scope for a PRNG algorithm is the safest solution to eschew lots of modifications tied to frequent hardware evolutions. This is the reason why we have decided to study in details the Saito proposition: MTGP.

DESCRIPTION OF MTGP

Data Structure

At the time we are writing this paper, MTGP was not referenced by any scientific publications yet, except that we found on the Internet that Saito's work is to be published soon [START_REF] Saito | A variant of mersenne twister suitable for graphics processors[END_REF]. We must describe its features in order to introduce the goal of our study.

First of all, MTGP is obviously inheriting from the properties of its elder, though it is quite different from a simple GPU implementation of the original MT, as seen in [START_REF] Podlozhnyuk | Parallel mersenne twister[END_REF]. As a matter of fact, Saito uses the original paper describing MT (Matsumoto and Nishimura, 1998) to lay his generator out. Thus, we will see that MTGP can suffer from some little problems already spotted for the MT "family". Since we often champion this family of generators, [START_REF] Reuillon | Simulations stochastiques en environnements distribués -Application aux grilles de calcul[END_REF] has studied 2 16 statuses of the original MT algorithm using the TestU01 Crush test battery from [START_REF] L'ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF]. The involved tests verify the linear complexity of the random sequence. MTGP is based upon the same linear recurrence to create random sequences, so it is not recommended for cryptographic purposes.

We have also noted that MTGP specified a common notion of the generators belonging to the parameterized family. We distinguish this cast of generators by the compound form of their data structures. It contains two distinct elements implied in the generation algorithm, we call them seed status and parameterized status. The first is basically the common seed given by the user to initialize a generator. The second stores parameters determined at a particular PRNG creation, it is supposed to be sort of a unique signature of the generator. Both these concepts were already present in MT. MTGP makes them more precise by explicitly using a data structure of the form we described in this paragraph. We propose a simple class diagram of this concept in figure 1: Parameterized statuses are a common way to ensure independence between parallel stochastic streams. This technique is called parameterization in the literature [START_REF] Hill | Practical distribution of random streams for stochastic high performance computing[END_REF]. Depending on the way it is settled, it can lead to poor results [START_REF] De Matteis | Controlling correlation in parallel monte carlo[END_REF]. Unfortunately, we do not have any mathematical theorem allowing us to check the independence between two generators, according to their parameters. However, MT came along with the DC algorithm, designed to create large sets of independent generators. This algorithm integrates an identifier, often the one of the processor or thread that will host the PRNG, to distinguish parallel random streams. It uses it as a part of the characteristic polynomial of the matrix used by the MT algorithm. We can easily conclude that if characteristic polynomials are prime to each other, their associated matrices will also be unique in the same context. Hence, we obtain independent parameter sets.

Hopefully this algorithm has been renewed for MTGP, enabling us to proceed in the same way. Furthermore, it improves the original algorithm by allowing the user to get a larger set of 2 32 parameters, where the original algorithm could only deliver 2 16 sets. This number is now too small for large computing grids such as the European Grid Initiative (EGI), with more than 240000 cores at the time we are writing this paper. This latter point forces us to keep skeptical concerning the independence of the MTGP produced by the new DC. Saito explains that a SHA1 (Secure Hash Algorithm) checksum of each characteristic polynomial is generated to let the user check he did not get duplicated entries.

Architecture Independence

One of the most interesting features of MTGP is to be available for both CPUs and GPUs architectures. Even if MTGP has been designed to run on GPUs, you can also find a CPU version at the same location (Saito, 2010a). We have based our study on the capability of the generator to merge transparently in CPU-based applications. Hence, we were able to test the PRNG on any CPU-based host with reliable and well-known tools like TestU01. This way, we have avoided the hazardous implementation of a new empirical test battery, which would have to be validated before.

Moreover, this property is really precious in our opinion. We intend to use this PRNG for stochastic simulations following the hybrid computing paradigm, where the sequential part of the application runs on a CPU host, while the parallel-one is executed on a GPU board. In such cases, stochastic streams will furnish random numbers to both the CPU and the GPU. With a PRNG like MTGP, we can keep our simulations homogenous, using the same PRNG on each computing element engaged. Handling independent parallel stochastic streams becomes understandably important when you have to deal with such hardware configurations. We will study in a next paragraph the specific DC coming along with MTGP to maximize this independence. Now considering the new elements we introduced in this whole section, we can extend our previous object model to the particular MTGP. Figure 2 depicts its main components: We have widely presented MTGP statuses in this section. Since this generator utilizes a PRNG per block of threads approach, we will need a different status per block to ensure the independence between the random streams produced. The sample program furnished by Saito takes a number of blocks to use as an argument and owns a set of 128 different statuses to feed these blocks. In the next section, we present the protocol that helped us to issue a large number of statuses for MTGP users.

MTGP BENCHMARK

Empiric Test of 10000 Statuses

We formerly introduced the DC tool, enabling us to create independent PRNGs to use in a parallel environment. It has supplied us 10000 independent parameter sets (parameterized statuses), each corresponding to a different MTGP. This step can be very computationally expensive when the algorithm has to look for generators displaying huge periods, such as 2 19937 for the original MT. According to (L'Ecuyer, 2010), periods contained between 2 100 and 2 200 should be sufficient for nowadays stochastic applications. Thus, we decided to manipulate generators of the lowest period allowed by the MTGP DC, which is still 2 3217 . Moreover, the lower the period is, the fewer it wastes bits to store the internal state vector of the PRNG, helping us to save some GPU memory. With this configuration, we have been able to get our entire ready to use MTGPs in a single day, using a 256-node Opteron-enabled Linux cluster.

The second phase consisted in applying the BigCrush test battery to each newly created generator, in order to check their quality. First of all, to easily analyze such an amount of results, we modified the TestU01 library output to enable it to produce lighter results output files. In this manner, we have been able to parse results files using script tools like Sed or Awk to generate statistics. Moreover, since lots of our computations have taken place on the European computing Grid Infrastructure (EGI), we reduced the quantity of data transferred on this slow bandwidth storage space. The use of this HPC tool was compulsory in our case, in fact [START_REF] L'ecuyer | TestU01:A Software Library in ANSI C for Empirical Testing of Random Number Generators -User's guide, detailed version[END_REF] forecasts BigCrush to take about 8 hours of CPU time on an average 64-bit processor. We could not afford to perform the equivalent of 80000 CPU-hours on a single cluster to get our results in a decent time.

Our final aim is to provide verified parameterized statuses to allow the simulation community to initialize GPU-enabled PRNGs without introducing any bias in stochastic simulations, according to the current knowledge. A basic selection consisted in keeping only statuses that had perfectly passed all tests of the battery. But given that this approach eliminated approximately 40% of the statuses, we tried to determine whether other statuses could be kept with a good confidence level. We set up a more formal analysis to answer this question.

Statistics-Based Analysis

Each test of the TestU01 Bigcrush battery [START_REF] L'ecuyer | Testu01: A c library for empirical testing of random number generators[END_REF] is governed by the H 0 hypothesis, that the successive output values of the RNG are i.i.d. U(0, 1), i.e. are independent random variables from the uniform distribution over the interval [0;1]. These tests are defined by a test statistic γ (which is a function of the numbers to be tested). They compute and report a number, called the p-value of the test, which is contained between 0 and 1. Furthermore, if γ has a continuous distribution, the p-value is U(0, 1) under H 0 . At this point, let us consider two precisions from L'Ecuyer and Simard:

1. « If the p-value is extremely small (e.g., less than 10 -10 , then it is clear that the RNG fails the test, whereas if it is not very close to 0 or 1, no problem is detected by this test. » (L'Ecuyer and Simard, 2007);

2. « Moreover, when a generator starts failing a test decisively, the p-value of the test usually converges to 0 or 1 exponentially fast as a function of the sample size when the sample size is increased further. » (L' Ecuyer and Simard, 2009).

According to these quotations, we decided to consider three p-value types, detailed hereafter:

• p-values contained between [0.001;0.999] are reckoned as correct, (these values are proposed by the TestU01 library);

• those included in the range [0;0,001[ U ]0.999;1] are counting as suspect;

• lastly, we refined the previous range since we needed to take into account extremely small p-values (less than 10 -10 ), called disastrous afterwards.

As mentioned previously, we ran our tests on 10000 independent statuses, with regards to MTGP DC. The column chart appearing on figure 3 represents the number of suspect p-values noticed for statuses where no disastrous p-values were obtained (due to layout considerations, only the thirty-two first tests are present on figure 3): The latter figure helped us to easily identify three test groups. They distinguish from others by recording more than 40 suspect p-values. The interesting point here is that these three groups characterize some aspects of the generator behavior. In fact, tests belonging to the same group are just differently parameterized versions of the same test. The noticeable tests are described as follows in the TestU01 documentation (L'Ecuyer and Simard, 2009):

• smarsa_CollisionOver (tests 9 to 12), is an overlapping pairs sparse occupancy (OPSO) test introduced in [START_REF] Marsaglia | A current view of random number generators[END_REF];

• snpair_ClosePairs (tests 22 to 25), is a m-nearest-pairs (m-NP) test (L'Ecuyer and Simard, 2009);

• swalk_RandomWalk1 (tests 74 to 79), applies simultaneously several tests based on a random walk of length l over the integers, for several (even) values of l (L' Ecuyer and Simard, 2009).

Under the H 0 ' hypothesis of a uniform distribution of the p-values over the interval [0;1], the distribution of the number of suspect p-values is binomial with parameters n = 10000 and p = 2/1000. So, we can reject the H 0 ' hypothesis with a very good confidence level (about 7.10 -6 for each test that accumulates more than 40 suspect p-values. However, we cannot do the same with H 0 . To do so, the γ-statistic used by the test should present a continuous distribution, which is not the case for the considered tests. Concretely, it means that pointing out suspect p-values brings useful information, but no matter the excesses of suspect p-values, such results do not provide a formal statistic proof of the test failure. That is why we use to consider a test is failed only when it returns disastrous p-values.

Six tests are missing on the previous figure: four, the 70, 71, 80, 81, were introduced in our MT description as problematic tests for any MT-like generator. So we increased the execution speed of the test battery simply by disabling those tests that would have systematically produced disastrous p-values. The other two tests, numbered 35 and 100, are the problematic ones. Only these assessments issue disastrous p-values in a non negligible quantity. We noted that about a tenth statuses were failing the 100 th Test, while more than a fifth went wrong with the 35 th Test. Figure 4 gives a graphical representation of the announced proportions: (Knuth, 1969). This test counts, for s = 0, 1, 2, . . . « the number of times that a sequence of exactly s successive values fall outside the interval [Alpha, Beta] (this is the number of gaps of length s between visits to [Alpha, Beta]). It then applies a chi-square test to compare the expected and observed number of observations. » (L' Ecuyer and Simard, 2009). A typical generator miscarrying this test « wanders in and out of [Alpha, Beta] for some time, then goes away from [Alpha, Beta] for a long while, and so on » (L' Ecuyer and Simard, 2007). However, one should note that the same test is perfectly passed with other Beta values. In view of the analysis we propose, this fact is obviously logical for Beta values lower than the incriminated one (1/32), but it is rather strange not to find disastrous p-values with a higher Beta value (Test 34 sets Beta to 1/16).

The test number 100 is referenced as sstring_HammingIndep test with N = 1, n = 107, r = 25, s = 5, L = 1200 and d = 0. It applies two tests of independence between the Hamming weights of successive blocks of L bits [START_REF] L'ecuyer | beware of linear congruential generators with multipliers of the form a = ±2q±2r[END_REF]. According to François Panneton, this test measures Hamming-weight dependencies between random values issued by a given generator. It tends to demonstrate that the recurrence does not shuffle bits enough from an iteration to another [START_REF] Panneton | Construction d'ensembles de points basée sur des récurrences linéaires dans un corps fini de caractéristique 2 pour la simulation Monte Carlo et l'intégration quasi-Monte Carlo[END_REF].

With this first experiment, we have put under the spotlight difficulties encountered by MTGP 2 3217 . Assuming that these problems are mostly concentrated on the two properties checked by tests 35 and 100, those producing disastrous p-values, we have focused our further studies on them.

PARAMETERIZED STATUS INFLUENCE

Seed Status Variation

Let us recall that the 10 000 tested statuses were only differing from their parameterized parts, whereas they shared the same seed status, fixed to 0. Viewing the previous results, we decided to work out whether the fact of passing a test or not, and by the way the quality of a generator, was due to either its parameterized status, or its seed status, or both. To do so, we settled a new experiment, sieving a set of 100 parameterized statuses alternately associated with 100 seed statuses randomly chosen. This initialization technique, called Random Spacing, represented a total of 10 000 combinations put to the proof of the guilty tests presented above. Figure 5 shows an extract of the graphical output for test 35. Red crosses mark a disastrous p-value, while blue ones indicate suspect p-values. An aggregation of red crosses on vertical lines shows that the parameterized status on the abscissa failed the considered test for all the seed statuses it was associated to. So, in the case of parameterized generators, it seems that the parameterized status establishes the generator quality on its own, independently from the selected seed status.

Period Variation

To comfort our previous assertion, we tried to make other elements of the parameterized status vary to observe their impact on the quality of the generator. As long as DC tries to figure out a tempering matrix such as the PRNG produces a well distributed sequence [START_REF] Matsumoto | Dynamic creation of pseudorandom number generators[END_REF], modifying this period should highly impact this property for the newly created statuses. Now, we previously brought forward that Test 35 (sknuth_gap) of BigCrush based its judgment on this characteristic. So, we intended to obtain much better results using 1000 MTGPs of period 2 23209 , confronted to tests 35 and 100 only. The result is crystal-clear since 99.5% of the statuses passed both tests without any problem. Moreover, we only noticed suspect p-values in the other 0.05%.

Obviously, a higher period eliminates sequence distribution issues, but this latter result could hide potential intrinsic weaknesses of the MTGP algorithm. Our last experiment will introduce as a standard a quality-proven PRNG of the same family: the original MT (Matsumoto and Nishimura, 1998).

Algorithm Variation

The original MT is designed with linear-recurrences preventing its recommendation for some particular applications such as cryptography. As far as we know, no other problems are referenced concerning this PRNG. That is why we consider it as a good standard to compare with the target of our investigations. An interesting point is that its own DC is able to produce parameters for the 2 3217 period, thus allowing us to work with MTs and MTGPs of the same period. This way, we focus our experiment on the algorithm-dependent parts of the parameterized statuses. Once again, we observed the behavior of each generator when faced to tests 35 and 100. We selected a sample of 1000 independent PRNGs to compare outputs with our previous benchmarks. Results are even more clear-cut than before: 99.9% of MT-dedicated statuses passed the two tests without any failure, whereas about 64% of MTGP statuses did.

The two previous results tend to show the influence of parameterized statuses. Here we have shown that this data structure is tightly bound to the algorithm using it. MTGP seems to present weaknesses when configured with shorter periods, since MT, running with the same relatively small period, eschews traps where its recent variant falls frequently into. This section results are summed up on the column chart displayed on figure 6:

CONCLUSION

This work was first intended to study a particular generator designed for GPU hardware architectures: MTGP. After an introduction to other recent approaches mentioned in the literature, in order to compare them with the MTGP strategy, we provided a description of MTGP and proposed a generic object model representing generators using distinct seeds and parameters. To complete our description, we achieved some analysis of this PRNG by facing it to the current most stringent test battery: BigCrush from the TestU01 test software library. Weaknesses identified during these experiments have been reduced by comparisons with other configurations of the generator as well as with the original MT.

The second goal of this study was to furnish parameterized statuses to the simulation community, to allow them to quickly and safely initialize their MTGPs without introducing any potential bias with a bad status. Only statuses that have passed BigCrush have been kept and will be freely available on the Internet.

We have shown that MTGP was safer with longer periods, according to TestU01 criteria, but the more a PRNG period is long, the more space it needs to store the internal state vector used by its algorithm. Nowadays, GPUs memory characteristics do not allow us to waste bytes to store PRNG data without influencing the whole application speed. By selecting only parameterized statuses referenced by our study or proofed by an equivalent benchmarking protocol, scientists using MTGP with the lowest available period (i.e. 2 3217 ) can dramatically reduce the memory footprint of their hybrid stochastic simulations. However, we should not forget that the empirical tests provided by TestU01 consider random sequences individually. So, we are now thinking of a way to test mixed random sequences, the way they can be in a single GPU application using several blocks of threads fed with different sources of randomness.

Figure 1 :

 1 Figure 1: Class diagram of a parameterized PRNG MTGP takes this idea one step further by introducing two kinds of statuses: the references and the fasts. The latters use pre-stored elements to decrease the initialization time, and programming techniques such as inlining to speed up the execution time. However, it results in a memory greedier status.

Figure 2 :

 2 Figure 2: Class Diagram for MTGP and its components

Figure 3 :

 3 Figure 3: Number of suspect results versus test numbers (extract displaying tests 1 through 32)

Figure 4 :

 4 Figure 4: Detailed Results for Tests 35 and 100 of the BigCrush battery

Figure 5 :

 5 Figure 5: Extracts of the results for Test 35: MTGP identifiers versus random seeds indices

Figure 6 :

 6 Figure 6: Percentage of passed results noticed for tests 35 and 100 depending on the PRNG

funded by the Auvergne Regional Council and the FEDER


Rebecca H Chisholm
email: rebecca.chisholm@inria.fr

Tommaso Lorenzi

Luis Almeida

Annette K Larsen

Alexander Lorz

Alexandre Escargueil

Jean Clairambault

Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, non-genetic instability and stress-induced adaptation

Keywords: Cancer modelling, reversible drug tolerance, phenotype evolution, non-genetic instability, stress-induced adaptation

In recent experiments on isogenetic cancer cell lines, it was observed that exposure to high doses of anti-cancer drugs can induce the emergence of a subpopulation of weaklyproliferative and drug-tolerant cells, that display markers associated with cancer stem cells. After a period of time, some of the surviving cells were observed to change their phenotype to resume normal proliferation, and eventually repopulate the sample. Furthermore, the drug-tolerant cells could be drug resensitized following drug washout. Here we propose a theoretical mechanism for the transient emergence of such drug tolerance.

In this framework, we formulate an individual-based model and an integro-differential equation model of reversible phenotypic evolution in a cell population exposed to cytotoxic drugs. The outcomes of both models suggest that selection, non-genetic instability, stress-induced adaptation and the interplay between these mechanisms can push an actively proliferating cell population to transition into a weakly-proliferative and drugtolerant state. Hence, the cell population experiences much less stress in the presence of the drugs and, in the long run, reacquires a proliferative phenotype, due to selection pressure and phenotypic fluctuations. These mechanisms can also reverse epigenetic drug tolerance following drug washout. Our study highlights how the transient appearance of the weakly-proliferative and drug-tolerant cells is related to the use of high-dose therapy. Furthermore, we show how stem-like characteristics can act to stabilize the transient, weakly-proliferative and drug-tolerant subpopulation for a longer time window. Finally, using our models as in silico laboratories, we propose new testable hypotheses that could help uncover general principles underlying the emergence of cancer drug tolerance.

Major findings

The application of our mathematical models to a recently studied experimental system (where there is evidence of an epigenetically regulated and transient drug-tolerant phenotype in an isogenetic cancer cell line), suggests that mechanisms based on fundamental laws of biology, namely a combination of selection, non-genetic instability mediated by phenotypic noise and stress-induced adaptation, can (reversibly) push an activelyproliferating and drug-sensitive cell population to transition into a weakly-proliferative and drug-tolerant state, which will eventually facilitate the emergence of more potent, proliferating and drug-tolerant cells.

Introduction

The interplay between intra-tumor heterogeneity [START_REF] Navin | Inferring tumor progression from genomic heterogeneity[END_REF][START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF] and perturbations in the tumor micro-environment [START_REF] Holohan | Cancer drug resistance: an evolving paradigm[END_REF][START_REF] Gupta | Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells[END_REF][START_REF] Song | Metformin kills and radiosensitizes cancer cells and preferentially kills cancer stem cells[END_REF][START_REF] Yang | Dynamic equilibrium between cancer stem cells and non-stem cancer cells in human sw620 and mcf-7 cancer cell populations[END_REF] plays a key role in the emergence of drug resistance in cancer cell populations. Intra-tumor heterogeneity results mainly from genetic modifications (e.g., mutations) leading, through a Darwinian-like process, to the selection of tumor cells expressing phenotypes adapted to the local environment [START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF][START_REF] Greaves | Clonal evolution in cancer[END_REF]. However, intra-tumor heterogeneity can also emerge from non-genetic processes mediated by either stochastic events or epigenetic modifications [START_REF] Glasspool | Epigenetics as a mechanism driving polygenic clinical drug resistance[END_REF][START_REF] Newman | Single-cell proteomic analysis of s. cerevisiae reveals the architecture of biological noise[END_REF][START_REF] Brock | Non-genetic heterogeneity-a mutation-independent driving force for the somatic evolution of tumours[END_REF]. In fact, perturbations of the tumor micro-environment induced by cytotoxic agents may 'instruct' cells to enter a more stress-resistant phenotype, which can be passed on to subsequent cell generations [START_REF] Pisco | Non-darwinian dynamics in therapy-induced cancer drug resistance[END_REF].

In recent experiments, performed on a genetically homogeneous population of lung (PC9) cancer cells, Sharma and co-workers [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF] showed that epigenetically regulated changes in phenotype can play an important role in the development of reversible drug tolerance. During these experiments, a small subpopulation of drug-tolerant cells was consistently detected, that could maintain viability in the presence of high-dose drug therapy.

These Drug-Tolerant Persisters (DTPs) were shown to be non-proliferative and display markers specific to Cancer Stem Cells (CSCs). After a period of time, approximately 20 percent of DTPs changed their phenotype to resume normal proliferation and lost CSC markers, still in the constant presence of the drugs. The resulting cells were labelled the Drug-Tolerant Expanded Persisters (DTEPs) (Fig. 1(a)). We may infer from these observations that two phenotypically distinct cell sub-clones subsequently emerge: a first one with a survival-prone but non-proliferative phenotype (the DTPs), and a second one, able to survive in extremely hostile conditions and proliferate (the DTEPs). Interestingly, both DTPs and DTEPs could be drug resensitized by drug-free passaging.

Therefore, the three distinct subpopulations-PC9s, DTPs and DTEPs-that compose the whole PC9 population at various times during drug treatment, although genetically identical, possess different functional phenotypes. Most notably, they can be characterized by their respective levels of survival potential in extreme environments and proliferation potential, which are two traits that have important clinical implications for cancer patients [START_REF] Aktipis | Life history trade-offs in cancer evolution[END_REF]. But what is driving the evolution of phenotype observed in the PC9 cancer cells? During drug therapy, the average level of survival potential (i.e., the average level of robustness towards life-threatening events in extreme conditions) in the PC9 cell population is non-decreasing over time. We argue that this can be explained by the interaction between two mechanisms. The first mechanism is non-genetic instability mediated by phenotypic noise [START_REF] Silander | A genomewide analysis of promoter-mediated phenotypic noise in escherichia coli[END_REF]. The term phenotypic noise, also known as biological noise [START_REF] Newman | Single-cell proteomic analysis of s. cerevisiae reveals the architecture of biological noise[END_REF][START_REF] Brock | Non-genetic heterogeneity-a mutation-independent driving force for the somatic evolution of tumours[END_REF], refers to the fluctuations in phenotype due to processes which may be intrinsic or extrinsic to the cell [START_REF] Swain | Intrinsic and extrinsic contributions to stochasticity in gene expression[END_REF], and may include, but are not limited to, noise related to gene expression [START_REF] Munsky | Using gene expression noise to understand gene regulation[END_REF][START_REF] Taniguchi | Quantifying e. coli proteome and transcriptome with single-molecule sensitivity in single cells[END_REF], protein production [START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF] and DNA methylation. Fluctuations in phenotype give rise to phenotype heterogeneity and can result in the transition of drug-sensitive cells into a robust phenotype [START_REF] Gupta | Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells[END_REF][START_REF] Brock | Non-genetic heterogeneity-a mutation-independent driving force for the somatic evolution of tumours[END_REF][START_REF] Mccullough | Plasticity of the neoplastic phenotype in vivo is regulated by epigenetic factors[END_REF]. These cells are able to survive drug exposure, and hence the second mechanism, selection, can cause the most robust cells to become the majority in the population.

The dynamics of the average level of proliferation in the cell population is even more interesting. Initially, as the majority phenotype transitions from the PC9 phenotype to the DTP phenotype, the average proliferation level drops. However, after a period of time, as DTEPs emerge from DTPs, the average level of proliferation increases back to normal.

The interplay between non-genetic instability and selection could explain the emergence of the DTEPs from the DTPs. Nevertheless, these mechanisms alone cannot account for the initial drop in the average level of proliferation of the surviving cell population. More specifically, if we assume there are no DTPs present in the initial PC9 population, then we would expect (small and rare) fluctuations in phenotype, combined with selection, to result in the emergence of DTEPs directly from PC9s (so that proliferation levels remain high). The transient dominance of DTPs indicates that another process is involved. A natural explanation is that there are a small number of DTPs present in the initial population of PC9 cells. These cells survive the initial drug therapy, and are subsequently selected for in the presence of the drugs. However, this scenario is only plausible if there are some DTPs present from the beginning. If this is not the case, then we propose that the DTP phenotype must be stress-induced [START_REF] Pisco | Non-darwinian dynamics in therapy-induced cancer drug resistance[END_REF][START_REF] Kreso | Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer[END_REF][START_REF] Williams | Quantitative high-throughput efficacy profiling of approved oncology drugs in inflammatory breast cancer models of acquired drug resistance and re-sensitization[END_REF]. Recent work on colorec-tal cancer cells by Kreso and co-workers [START_REF] Kreso | Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer[END_REF] showed that, in response to chemotherapy, actively proliferating cells were preferentially eliminated compared with other relatively dormant progeny. The authors suggest that these cancer cells may use dormancy (i.e., mere survival without differentiation, senescence nor commitment into the division cycle) as an adaptive strategy. Accordingly, we suggest that, if there are no DTPs present in the initial PC9 population, then stress induces the PC9 cells to adopt a more DTP-like phenotype [START_REF] Pisco | Non-darwinian dynamics in therapy-induced cancer drug resistance[END_REF][START_REF] Kreso | Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer[END_REF]. In this adaptive scenario, PC9 cells are induced to redistribute their finite energy resources away from proliferation-oriented tasks and towards survival-oriented tasks (for example, nucleotide metabolism, DNA synthesis and repair [START_REF] Lackner | Mechanisms of acquired resistance to targeted cancer therapies[END_REF]). On the other hand, robust cells, with a high survival potential, experience much less stress in the presence of the drugs. They have no need to redistribute their energy resources away from proliferation-oriented tasks. Accordingly, in response to drug exposure, one will expect a population of mostly drug-sensitive PC9 cells to experience an initial sharp drop in average proliferation levels. Only when levels of resistance towards drug damage increase sufficiently, due to fluctuations in phenotype and selection, will cells be able to reallocate resources back into proliferation-oriented tasks, and resume normal proliferation [START_REF] Davidson | An increase in the expression of ribonucleotide reductase large subunit 1 is associated with gemcitabine resistance in non-small cell lung cancer cell lines[END_REF][START_REF] Gautam | Suppression of lung tumor formation by the regulatory subunit of ribonucleotide reductase[END_REF].

Motivated by these considerations, here we propose an Individual-Based (I-B) computational model [START_REF] Pisco | Non-darwinian dynamics in therapy-induced cancer drug resistance[END_REF][START_REF] Anderson | Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment[END_REF][START_REF] Charlebois | Gene expression noise facilitates adaptation and drug resistance independently of mutation[END_REF][START_REF] Silva | Research a theoretical quantitative model for evolution of cancer chemotherapy resistance[END_REF] and an Integro-Differential Equation (IDE) model [START_REF] Delitala | A mathematical model for the dynamics of cancer hepatocytes under therapeutic actions[END_REF][START_REF] Lavi | The role of cell density and intratumoral heterogeneity in multidrug resistance[END_REF][START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF] of the phenotype evolution observed in the cell lines of Sharma and co-workers [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF]. Such models can be used as in-silico laboratories to highlight stylized facts, and uncover mechanisms that underlie emergent features of cancer cell populations. The I-B computational model allows an intuitive and flexible description of the system at hand, while the IDE model makes it possible to study the system in terms of qualitative and asymptotic analysis, and is computationally less expensive.

In particular, we focus on the evolution of two phenotypic traits which show substantial variability during drug treatment and after drug washout in the PC9 cell lines: a cell's survival potential and proliferation potential [START_REF] Aktipis | Life history trade-offs in cancer evolution[END_REF]. These models rely on the assumption that PC9 cells possess significant phenotypic plasticity. Furthermore, we assume that, during drug treatment, the surviving cell population can consist of cells residing within a spectrum of intermediate phenotypic states, ranging from PC9 through to DTP, and to DTEP. This idea is consistent with observations in other genetically identical cancer cell populations, where there exists significant functional phenotype heterogeneity, even in homogeneous environments [START_REF] Gupta | Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells[END_REF][START_REF] Pisco | Non-darwinian dynamics in therapy-induced cancer drug resistance[END_REF][START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF][START_REF] Kreso | Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer[END_REF]. For example, although epithelial-mesenchymal transition is often portrayed as a switch between two distinct states, in fact it has been shown that epithelial cells progressively acquire mesenchymal features at the same time as they lose epithelial attributes [START_REF] Chao | Partial mesenchymal to epithelial reverting transition in breast and prostate cancer metastases[END_REF]. In the same way, during drug treatment, we assume that the PC9s can gradually acquire drug tolerance at the same time as they gradually change their proliferation ability.

In the presence of high-dose drug therapy, we illustrate the ability of the models to recover the evolution of the PC9 cells into DTPs and finally into DTEPs, and also the drug resensitization of DTPs and DTEPs. Showing this, we can then argue that the models have validity for suggesting plausible evolutionary mechanisms underlying drug tolerance. Therefore, we utilize our I-B computational and IDE models in parallel to address the following questions: Q1. Is non-genetic instability necessary in the development of drug tolerance? Q2. How do the evolutionary dynamics of the cell population compare between the regimes of low and high drug dose? Q3. Is it possible to determine whether therapy induces the DTP phenotype? Q4. What is the role of stem-like characteristics of the DTP cells in the development of drug tolerance?

Our study highlights the important role of non-genetic instability in the emergence of DTEPs, and also leads us to propose experiments which may help to assess whether there are DTPs present ab initio, and whether PC9 cells are undergoing stress-induced adaptation of their phenotype. Finally, our models predict that the transient dominance of the DTP subpopulation is regulated by their stem-like characteristics and is a direct result of the high doses of chemotherapeutic drugs used in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF].

Quick guide to equations and assumptions

We describe the evolution of phenotype in a well-mixed PC9 cancer cell population exposed to cytotoxic drugs (whose concentration at time t is c(t)), using both an I-B and an IDE formalism. We consider only the case of constant infusion of cytotoxic drugs (i.e., c(•) := c > 0) and characterize the state of each PC9 cell by its level of expression of two phenotypic traits: survival potential and proliferation potential. In this framework, we identify a PC9 cell as having a low value of survival potential and a high value of proliferation potential, a DTP cell as having a high value of survival potential and a low value of proliferation potential, and a DTEP cell as having a high value of survival and proliferation potential (schematized as in Fig. 1(b)). We model the three key biological mechanisms of phenotype evolution in the PC9 cancer cells, namely selection, stress-induced adaptation and non-genetic phenotype instability, by using three separate mathematical strategies, which are tailored to fit the I-B and the IDE formalisms. Specifically, we model natural selection through a proliferation probability p and a death probability d, which depend on the levels of proliferative and survival potentials of the cells, as well as on the cell micro-environment at time t. Stress-induced adaptation of cell proliferation level, which occurs in the presence of the drugs, and depends on the level of survival potential, is modeled as an advection term. It leads to a decrease in the level of proliferation, and varies with the survival potential of the cell and the concentration of the drugs. Finally, non-genetic phenotype instability is modeled by stochastic fluctuations, or diffusion, in the levels of cell proliferation and survival potentials.

Individual-based model

In the I-B formalism, we consider each cell as an individual agent, and label it by an index 0 ≤ i ≤ N (t), where N (t) ≥ 0 is the size of the population at time t ∈ [0, T ], and T is the end time of the simulation. The normalized expression levels of the survival-potential and proliferation-potential traits, in each cell i at time t, are modeled, respectively, by the random variables

X i (t) : [0, ∞) → [0, 1] and Y i (t) : [0, ∞) → [0, 1].
We simulate the evolution of the cell population in discrete time, according to the algorithm pictured in Supplementary Fig. S2. Over the time interval between two successive time instants t and t + ∆t, we first allow each cell i either to proliferate, undergo apoptosis or remain in a quiescent state according to the respective probabilities p(X i , Y i , N (t))∆t, d(X i , c(t))∆t and 1 -∆t (p(X i , Y i , N (t)) + d(X i , c(t))). If a cell proliferates, we assume both daughter cells inherit the parent's trait values. This is consistent with experimental observations in cancer cell lines, which suggest that a cell's "epigenetic memory" can last multiple cell generations [START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF]. After all cells have undergone one iteration of the proliferation and death process, and returned to the quiescent state, we then let each cell update its trait values according to the following system of discretized Stochastic Differential Equations:

X i (t+∆t) = X i (t)+D √ ∆t W 1 i (t), non-genetic instability Y i (t+∆t) = Y i (t)+D √ ∆t W 2 i (t) non-genetic instability + ∆t v(X i (t), c(t); v),
stress-induced adaptation of the proliferation level [START_REF] Navin | Inferring tumor progression from genomic heterogeneity[END_REF] where v represents the average sensitivity of the cell proliferative potential to stressinducing agents, and D is the average rate of phenotypic fluctuations.

Integro-differential equation model

In the IDE formalism, we consider the cell population to be structured by two continu-

Extended individual-based model

In the extended I-B model, we keep the birth and death process unchanged from the original I-B model. However, instead of allowing the DTPs to update at each time step, we consider them to be stem-like so that they do not change their phenotype, unless they proliferate. When a DTP undergoes proliferation, we allow its offspring to update according to p 1 , p 2 and p 3 which represent, respectively, the probabilities for symmetric self-renewal, asymmetric self-renewal and symmetric differentiation. The update mechanisms of the PC9s and DTEPs remain unchanged. This new algorithm is pictured in Supplementary Fig. S3.

The assumptions on the functions p, d and v are summarized in Supplementary Tables S1-S2. All definitions of the parameter functions are provided in the Supplementary Material.

Results

For each scenario presented below, we consider the concentration of cytotoxic drugs (i.e., the parameter c) to be high if it is greater than, or equal to, the LC 90 dose (defined as the drug dose required to kill 90% of the total cell population in the initial stage of drug therapy, before the population starts to recover). On the other hand, we assume that a drug dose is low if it is smaller than, or equal to, the LC 50 dose (defined as the drug dose required to kill 50% of the total cell population in the initial stage of drug therapy, before the population starts to recover). The end-time of the simulations is represented by the parameter T , and is defined as the time required for the global population density to reach equilibrium.

Mathematical models recover the population level dynamics of the PC9 cells

To determine whether our models can reproduce the experimental observations of the PC9 cancer cell lines during drug therapy, as reported in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], we numerically solve the IDE model and simulate the I-B model for a high dose of the cytotoxic drugs. Both models can successfully reproduce the important population level dynamics of the PC9 cells, namely the extinction of the PC9 cells, followed by the transient appearance of DTPs and the emergence of the DTEP subpopulation.

In Fig. 2, we present numerical solutions from the I-B and IDE models for the case when there are a few DTPs present in the initial population (i.e., 5% of the global population density), and when there is no stress-induced adaptation of the proliferation potential. In the experiments reported in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], both the resistant DTPs and DTEPs can be drug resensitized by drug-free passaging. To test whether this is also the case in our models, we numerically solve the IDE model and simulate the I-B model for a period of time with high doses of cytotoxic drugs, so that either DTPs or DTEPs are left as the dominant subpopulation. We subsequently switch the drug concentration off to simulate instability is present, the long-term total global population density can be much smaller.

In fact, if the concentration of the drugs is high enough, then the population can go extinct when non-genetic instability is absent. These predictions are illustrated in Fig. 4(c)-(d) (see also Supplementary Fig. S5(c)-(d)).

Therefore, we hypothesize that non-genetic phenotype instability is essential for the emergence of DTEPs, and for the development of drug tolerance in PC9 cells. Furthermore, the absence of non-genetic instability can result in the stabilisation of the DTP phenotype in the surviving population, and even extinction of PC9 cells during high-dose drug therapy.

A2. A high dose of cytotoxic drugs is necessary for the transient dominance of

DTPs

We use our IDE model to investigate the evolutionary dynamics of the PC9 cells in a low drug-dose regime. For Scenario (i), described above, our model predicts that, ceteris paribus, decreasing the drug concentration can prevent the emergence of DTPs and speed up population recovery. These dynamics are illustrated in Fig. 5(a)-(b), where we plot the trajectories of the mean phenotypic expression of the population and the corresponding global population density as a function of time, for both a high dose and a low dose of the cytotoxic drugs. This is also true for Scenario (ii) (see Fig. 5(c)-(d)). Therefore, we can conclude that a higher dose of the cytotoxic drugs is a key ingredient to see the transient dominance of DTPs.

A complete set of comparisons between the high and low drug-dose regimes, for different values of v (the average sensitivity of the cell proliferative potential to stress-inducing agents) and β (the average rate of phenotypic fluctuations), are provided in Figures S5-S6

and support this conclusion.

A3. Low cytotoxic drug-dose experiments could establish whether the DTP phenotype is stress-induced

For the case when there are a few DTPs present in the initial PC9 population, our model predicts that in a low drug-dose regime, the emergence of a meta-stable subpopulation of DTP cells can require the presence of stress-induced adaptation, non-genetic instability in phenotype, and the proper interplay between them. This is illustrated in Fig. 6(b) where, for a given value of β (the average rate of phenotypic fluctuations) the trajectories of the population mean trait levels pass through the region corresponding to the DTP phenotype, only when stress-induced adaptation is present. In fact, in order to allow for the emergence of a metastable subpopulation of DTP cells, a sufficiently high sensitivity of the cell proliferation potential to stress-inducing agents (i.e., the parameter v) and a suitable value of β are required. This is in contrast to the high drug-dose regime where DTPs emerge in the absence of stress-induced adaptation, and the evolutionary dynamics of the population are unaltered by changes in v (see Supplementary Fig. S6(c)-(d)).

Another interesting prediction from our model in the low-dose regime is that, if there are a few DTPs present in the initial PC9 population, and stress-induced adaptation is absent, then it is possible to see an increase in the average proliferative potential in the surviving population, followed by a decrease before a second increase (see Fig. 6). This is not predicted in the high-dose regime, nor is it predicted if there are no DTPs present initially, and it is related to the timing of the emergence of the DTEPs, and the disappearance of the PC9s.

Therefore, one way to determine whether there are pre-existing DTPs within the initial PC9 population, and whether stress-induced adaptation is occurring, would be to repeat the experiments of [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF] with a low dose of the cytotoxic drugs, and observe the average mitotic rate of the surviving population. Our models predict that, if there is an initial increase in the average mitotic rate, followed by a decrease and a further increase, then it is likely that DTPs are present in the initial population, and stress-induced adaptation is not occurring.

A4. Stem-like characteristics temporarily stabilize the subpopulation of DTPs

Sharma and co-workers showed that the DTPs display stem-like surface markers, which are absent in the DTEPs [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF]. Therefore, it is reasonable to assume that the DTPs express some stem-like characteristics. In particular, stem cells have the potential to generate more stem cells, as well as differentiated daughter cells [START_REF] Morrison | Asymmetric and symmetric stem-cell divisions in development and cancer[END_REF]-they are capable of symmetric cell divisions, where a dividing stem cell can produce either only stem-cell daughters (selfrenewal) or only differentiated daughters, or asymmetric cell divisions, where a dividing stem cell produces one stem-cell daughter and one differentiated daughter. Under steadystate conditions, it is generally thought that stem cells divide asymmetrically. However, when a population of stem cells is subjected to chemotherapy, they have the ability to switch to a symmetric mode of division (symmetric self-renewal), in order to regenerate depleted stem-cell pools [START_REF] Morrison | Asymmetric and symmetric stem-cell divisions in development and cancer[END_REF][START_REF] Bodine | Bone marrow collected 14 days after in vivo administration of granulocyte colony-stimulating factor and stem cell factor to mice has 10-fold more repopulating ability than untreated bone marrow[END_REF].

Motivated by these considerations, we utilize the extended I-B model to investigate the role of self-renewal and asymmetric cell divisions in the development of drug tolerance.

We analyse the dynamics of the cells for a range of different values of the symmetric self-renewal probability p 1 , when the probability of asymmetric self-renewal p 2 equates to 1 -p 1 , so that the probability of symmetric differentiation p 3 is zero (see Supplementary Fig. S3).

Again we consider the two hypothetical scenarios described in Section A1. In both scenarios, if DTPs do transiently dominate the surviving population, then a higher probability of self-renewal increases the time taken for the emergence of the DTEPs from the DTPs. Furthermore, after the DTEPs have emerged, a higher probability of self-renewal corresponds to a larger and more stable subpopulation of DTPs in the surviving population. These predictions are illustrated in Fig. 7 where, for Scenario (i), we plot the trajectories of the mean phenotypic expression of the population and the corresponding numbers of PC9s, DTPs an DTEPs as a function of time, for three different values of p 1 (see also Supplementary Fig. S8(c)-(d) for Scenario (ii)).

On the other hand, if the DTPs do not transiently dominate the surviving population, either because the population becomes extinct before it can become drug tolerant (see Supplementary Fig. S7(c)-(d)), or because the DTEPs emerge directly from the PC9s (see Fig. S8(a)-(b)), then, not surprisingly, the probability of self-renewal p 1 does not play any role in the dynamics of the population.

Discussion

The PC9 cell line experiments reported in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF] are performed in an isolated and relatively homogeneous environment, and involve only a few constituents-a genetically identical cell population, culture media, and combinations of drugs. Furthermore, each experiment has clear observables, namely the percentage of surviving cells and their phenotypic distributions, which make the PC9 cancer cell lines an ideal system to study from a mathematical perspective.

Here, we have presented an I-B computational model and an IDE model of the evolution of phenotype observed in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], which rely on the assumption that this evolutionary process is driven by the interplay between selection, non-genetic phenotype instability and stress-induced adaptation. Distinct mathematical strategies have been utilized to model each evolutionary mechanism. Specifically, we have modeled selection through prolifer-ation and death probabilities, which depend on the levels of proliferative and survival potentials of the cells, as well as on the cell micro-environment. Non-genetic phenotype instability has been modeled by stochastic fluctuations, or diffusion, in the levels of cell proliferation and survival potentials. Finally, stress-induced adaptation of cell proliferation level, which occurs in the presence of the drugs, and depends on the level of survival potential, has been modeled with an advection term. This enabled us to better understand the effect of each mechanism on the phenotypic evolution in the PC9 cell line.

The mathematical models presented here reproduce the main experimental observations detailed in [START_REF] Swain | Intrinsic and extrinsic contributions to stochasticity in gene expression[END_REF], and also provide insight into the mechanisms controlling the dynamic variation in phenotype of the PC9 cancer cells during drug therapy, and after drug washout. Using our models as in silico laboratories, we performed mathematical analyses with two types of initial conditions, one where the population consists of only PC9s, and the other where there are also a small number of DTPs. Our analysis highlights the important role of non-genetic fluctuations in phenotype in the emergence of drug tolerance in PC9 cancer cell lines. In particular, we suggest that the absence of non-genetic instability can result in the stabilisation of the DTP phenotype in the surviving population, so that DTEPs do not emerge, or even in extinction. This is a key result since it supports the idea that epigenetic therapy may be a promising therapeutic strategy in the war against cancer [START_REF] Clozel | Mechanism-based epigenetic chemosensitization therapy of diffuse large b-cell lymphoma[END_REF][START_REF] Juergens | Combination epigenetic therapy has efficacy in patients with refractory advanced non-small cell lung cancer[END_REF][START_REF] Timp | Cancer as a dysregulated epigenome allowing cellular growth advantage at the expense of the host[END_REF]. Furthermore, it is also consistent with the experiments reported in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], where epigenetic drugs, when used in combination with high doses of the cytotoxic therapy, were able to prevent the emergence of drug tolerance and kill the cancer cell population.

Another important prediction of our models is that the transient dominance of DTPs is strictly related to the usage of high doses of cytotoxic drugs. If experimentalists apply a lower dose of cytotoxic agents to the PC9 cell population during drug therapy, we propose that it would be highly unlikely to observe DTPs. Rather, we would expect the DTEPs to emerge directly from the PC9 population. Note that this is the usual way to yield stable drug-tolerant lineages [START_REF] Gottesman | Mechanisms of cancer drug resistance[END_REF].

Our analyses also suggest that stem-like, self-renewing cell divisions of DTPs can act to increase the number of DTPs in the system, and is consistent with other mathematical models of stem-cell dynamics [START_REF] Marciniak-Czochra | Modeling of asymmetric cell division in hematopoietic stem cells-regulation of self-renewal is essential for efficient repopulation[END_REF]. In particular, we expect that DTP self-renewal can act to stabilize the DTP subpopulation for a longer period of time during drug therapy, compared to the case without self-renewal, as well as to delay the emergence of DTEPs.

Based on our models, we can conclude that if there are no DTPs present in the initial population of PC9 cells, then it is likely that a proper interplay between selection, non-genetic phenotype instability and stress-induced adaptation is mandatory for the transient appearance of the DTP phenotype during high-dose drug therapy. On the other hand, if there are some DTPs present in the initial population, and the concentration of the drugs is high enough, then selection and phenotypic fluctuations are enough to explain the experimental observations reported in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF]. Therefore, the next biologically meaningful question is: are DTPs present in the initial population? Our analysis enabled us to propose a low cytotoxic drug-dose experiment which could answer this question.

Another possible way to determine whether DTPs are present ab initio would be to utilize the experimental technique described in [START_REF] Moore | Slow-cycling therapy-resistant cancer cells[END_REF] to identify and isolate slow-cycling cells in the parental PC9 cell population, and subsequently compare the drug sensitivity of these cells with that of the more rapidly dividing PC9 cells. The presence of slowcycling, drug-tolerant cells would suggest that DTPs are present in the initial population of PC9s.

The low-dose experiment we proposed could also determine whether stress-induces the DTP phenotype. However, another way to assess whether PC9s are induced by stress to acquire the DTP phenotype would be to perform a fluorescence-activated cell sorting analysis on the PC9s with respect to the Cluster of Differentiation 24 (CD24) surface marker, at a time soon after drug delivery. In [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], the authors showed that the DTPs have a high level of expression of CD24, compared with the PC9s and DTEPs.

Therefore, if the surviving population has an intermediate level of expression soon after drug exposure, relative to the expression of CD24 in PC9s and DTPs, then this would suggest that the DTP phenotype is a result of non-genetic fluctuations in phenotype, selection and stress-induced adaptation.

Finally, we believe that structuring cancer cell populations according to a continuous phenotype is, in general, an appropriate strategy to describe the evolutionary dynamics of cells under stress. This is particularly true for the cell population studied in [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], where phenotypes are a clear source of biological variability. In our study, we have considered a two-dimensional structuring phenotype. However, we could introduce a third phenotypic trait to represent (de-)differentiation, which may help to better characterize the evolution of the system, particularly for what concerns the dynamics following drug washout.

It is then natural to wonder about what could be a biological quantitative variable underlying such a continuous (multidimensional) phenotype. We propose, as a reasonable candidate, the degree of (transient) epigenetic modifications of the DNA (methylation, histone acetylation) in relevant parts of the genome. Such modeling has been previously

proposed in [START_REF] Lei | Mathematical model of adult stem cell regeneration with cross-talk between genetic and epigenetic regulation[END_REF], using a one-dimensional, phenotype-like structure variable. By using a multidimensional phenotype, we propose to consider differential effects of epimutations on genes responsible for the classical intracellular pathways of proliferation, survival and (if necessary, in a three-dimensional extension of the present model) differentiation.

( 2 )

 2 where β is related to the parameter D from the I-B model through the relation β = D 2 /2.

Fig. 2 (

 2 Fig. 2(a) displays the trajectories of the mean phenotypic expression of the population.It is clear that, as time progresses, the mean phenotypic expression changes from being PC9-like (high proliferation and low survival potentials), to being DTP-like (low proliferation and high survival potentials) and finally to being DTEP-like (high proliferation and survival potentials). Such dynamics are also shown in Fig.2(c)-(j), which displays the phenotypic distribution at different times during drug therapy. In Fig.2(b), the plot of the global population density as a function of time reveals that, shortly after drug exposure, there is a large reduction in cell number (more than 90%), followed by a time of relatively constant population level, before the population recovers back to normal levels.Similar results can be obtained when there are no DTPs present in the initial population, and when we include stress-induced adaptation of proliferation levels (we present numerical solutions for these cases, from the IDE model in Fig.4(a)-(b) and the top rows of Figures S5-S6).

Fig. 2

 2 Fig.2also shows that we can obtain a good match for the dynamics predicted by the I-B model with that predicted by the IDE model. We also achieve the same qualitative match for the other cases considered in this section. Therefore, with the exception of our study on the stem-cell characteristics of DTPs (where we use the extended I-B model), from now on we present numerical results obtained from the IDE model only.

Figure 1 .

 1 Figure 1. Phenotype evolution in PC9 cells during cytotoxic-drug therapy. (a) Schematic diagram of phenotype evolution in PC9 cells during high-dose cytotoxic-drug therapy. Note that, in the scenario pictured here, there are DTPs present in the initial PC9 population. However, this may not necessarily be the case. (b) Schematic diagram representing the evolution of proliferation and survival potential levels in PC9 cells during cytotoxic-drug therapy. The dashed-red lines highlight the regions of the phenotype space corresponding to the PC9s, DTPs and DTEPs.

Figure 2 .

 2 Figure 2. Mathematical models recover the population level dynamics of the PC9 cells. (a) Trajectories of the mean phenotypic expression for 0 ≤ t ≤ T . (b) Global population density as a function of time. (c) Cell distribution over the survival potential at three different times t = T /50, T /10, T /2. (d) Cell distribution over the proliferation potential at three different times t = T/50,T/10,T/2. Solid lines are the numerical results from the IDE model, while the broken lines are from the I-B model. (e)-(g) Phenotypic distribution of the cell population at three different times t = T /10, T /2, T for the I-B model. (h)-(j) Phenotypic distribution of the cell population at three different times t = T /10, T /2, T for the IDE model. Details are provided in the Supplementary Material.

Figure 3 .

 3 Figure 3. Mathematical models recover the drug resensitization of (a) DTPs and (b) DTEPs following drug washout. Trajectories of the mean phenotypic expression. Here, solid lines are the numerical results from the IDE model during drug exposure, while the broken lines are from the IDE model after drug washout. Details are provided in the Supplementary Material.

Figure 4 .

 4 Figure 4. Non-genetic instability is crucial for the emergence of DTEPs. We illustrate the effect of the rate of phenotypic fluctuations (i.e., the parameter β) on (a), (c) the trajectory of the population mean trait levels for 0 ≤ t ≤ T , and (b), (d) the corresponding global population density, as a function of time, in a high drug-dose regime and when there are (a), (b) only PC9s present in the initial population and stress-induced adaptation is present, or (c), (d) 98% PC9s and 2% DTPs present in the initial population and stress-induced adaptation is absent. Note that in panels (a) and (c), the trajectories corresponding to the case β = 0, terminate (when the populations become extinct) in the regions corresponding, respectively, to the PC9 and DTP phenotype. Details are provided in the Supplementary Material.

Figure 5 .

 5 Figure 5. A high dose of cytotoxic drugs is necessary for the transient dominance of DTPs. We illustrate the effect of the concentration of the cytotoxic drugs (i.e., the parameter c) on (a), (c) the trajectory of the population mean trait levels for 0 ≤ t ≤ T (where T is the end time of the simulation), and (b), (d) the corresponding global population density, as a function of time, when there are (a), (b) only PC9s present in the initial population and stress-induced adaptation is present, or (c), (d) 98% PC9s and 2% DTPs present in the initial population and stress-induced adaptation is absent. Details are provided in the Supplementary Material.

Figure 6 .

 6 Figure 6. Low drug-dose experiments could establish whether the DTP phenotype is stress-induced. Trajectories of the population mean trait levels during low drug-dose therapy, when there are 98% PC9s and 2% DTPs present in the initial population, and

Figure 7 .

 7 Figure 7. Stem-like characteristics of DTPs stabilize the transiently dominant subpopulation of DTPs for longer. We illustrate the effect of the DTP self-renewal probability on (a) the trajectories of the mean phenotypic expression and (b) the corresponding number of PC9s (solid lines), DTPs (dotted lines) and DTEPs (dash-dot lines) as a function of time. Numerical results come from the extended I-B model. Details are provided in the Supplementary Material.

Figures

  Figures

Figure 1

 1 Figure 1

Figure

  Figure 3

3

 3 Figure 3

Figure 4 Figure 5

 45 Figure 4

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

  ous, real variables x ∈ [0, 1] and y ∈ [0, 1], which represent, respectively, the normalized expression levels of the survival-potential and proliferation-potential traits. The popula-

				1 0	1 0 n(x, y, t) dxdy, and the
	time evolution of n is governed by the following equation:
	∂n ∂t	(x, y, t) +	∂ ∂y	v(x, c(t); v)n(x, y, t)
				stress-induced adaptation
				of the proliferation level

tion density of cancer cells is modeled by the function n(x, y, t) ≥ 0, where the global population density at time t ∈ [0, T ] is computed as (t) = = p(x, y, (t)) -d(x, c(t)) n(x, y, t) selection + β ∆n(x, y, t), non-genetic instability

Acknowledgments

The research leading to this paper was (partially) funded by the french "ANR blanche" project Kibord: ANR-13-BS01-0004. TL was supported by the Fondation Sciences Mathématiques de Paris (FSMP) and by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference: ANR-10-LABX-0098).

drug washout, and allow the system to continue to evolve. We find that the models can successfully reproduce drug resensitization of the DTPs and DTEPs (see Fig. 3). This result also holds after drug washout following low-dose drug therapy (see Supplementary Fig. S4).

A1. Non-genetic instability is crucial for the emergence of DTEPs

To investigate the role of non-genetic instability in the evolution of drug tolerance in PC9 cells, we analyse the phenotype evolution in the PC9 population, for different values of the average rate of phenotypic fluctuations (i.e., the parameter β) for a high dose of cytotoxic drugs. We consider two hypothetical scenarios: Scenario (i) when the initial population consists of only PC9s and stress-induced adaptation of the proliferation level is present; Scenario (ii) when a small number of DTPs are present initially and stress-induced adaptation of the proliferation level is absent.

For the first scenario, our model predicts that if the average rate of phenotypic fluctuations is too large, then the PC9s will evolve directly into the DTEPs. Furthermore, if this rate is too small, then the PC9s will become extinct before they are able to express the more robust DTP phenotype. These predictions are illustrated in Fig. 4(a)-(b), where we plot the trajectories of the mean phenotypic expression of the population and the corresponding total cell density as a function of time, for three values of β (see also Supplementary Fig. S5(a)-(b)).

For the second scenario, the model predicts that in the absence of non-genetic instability, the DTEPs cannot emerge from the DTPs. Therefore the DTP phenotype is stabilized inside the surviving population. Consequently, compared with the case when non-genetic