José Sanchez
email: jose.sanchez@g2elab.grenoble-inp

Raphaël Caire
email: raphael.caire@g2elab.grenoble-inp

Nouredine Hadjsaid
email: nouredine.hadjsaid@g2elab.grenoble-inp

Towards a Complex Networks Modeling of Interdependent Critical Infrastructures

La fiabilité et la sécurité des Infrastructures Critiques mettent l'accent sur la nécessité d'identifier et de comprendre les vulnérabilités qui émergent de l'interaction entre eux, notamment entre les systèmes TIC et les réseaux électriques. La modélisation des relations physiques, logiques et géographiques peuvent donner quelques informations sur les interactions entre les infrastructures pour éclairer leurs interdépendances et identifier leurs vulnérabilités. Cet article présente une méthode, developé dans le cadre du Projet SINARI, qui adapte la théorie des nombres complexes à la théorie des réseaux complexes. Le résultat de cette symbiose est un modèle qui permet d'identifier les vulnérabilités inhérentes des infrastructures couplés. La méthode est démontrée sur un réseau de distribution typique français, y compris un système TIC.

Introduction

Systems are all around people, and each day the connections between them are getting increasingly complex. In fact, the world even depends on some of these systems for its survival [START_REF] Rinaldi | Indentifying, understanding, and analyzing critical infrastructure interdependencies[END_REF]. Governments use the term "Critical Infrastructures" to refer to systems that are essential to the defense and economic security of their nations. The Commission of the European Communities identified 9 key Critical Infrastructures, including Energy networks, Information and Communication Technologies (ICT) and government national sites and monuments [START_REF]Critical Infrastructure Protection in the Fight against terrorism[END_REF].

One important infrastructure concerns the Energy Network, which includes the production, refining, storage, and distribution of oil, gas, and electric power. However, this infrastructure integrates high-speed, reliable, bidirectional, and secure data communications networks. Particularly, Distribution Domain integrates external communication flows with control centers, electricity markets, transmission network, end-users, distributed storage, and distributed generation; and internal communication flows, i.e. between control, measure, and protection components. These interactions are presented in Fig. 1 [START_REF]NIST Framework and roadmap for smart grid interoperability standards[END_REF]. This fusion of Power Grid with ICT has added complexity to an already complex field. Recent events [START_REF] Falliere | W32. Stuxnet dossier[END_REF]- [START_REF] Matrosov | Stuxnet under the microscope (v1.0)[END_REF] have shown that failures in one infrastructure can affect other infrastructures and that it is not possible to protect an infrastructure without identifying and understanding its vulnerabilities, i.e. multi-infrastructure threats. As a result, when talking about vulnerability studies in the Power Systems, one cannot consider a single system, but a system connected to other systems.

FIG. 1 : Distribution Domain (NIST Smart Grid

Framework 1.0) Although various researchers have studied multiinfrastructures vulnerability, those studies are still in an early stage and many questions remain unanswered [START_REF] Kröger | Critical infrastructures at risk: a need for a new conceptual approach and extended analytical tools[END_REF]. For example, Petri Nets models [START_REF] Tranchita | ICT and power systems: An integrated SINARI page 6/6 approach[END_REF] showed that, despite the stochastic analysis they perform, these models are impractical for large systems due to the great amount of manual effort and the lack of flexibility to evaluate different scenarios. Some of the unanswered questions also revolve around the mutual behavior of connected infrastructures; as most proposed methods evaluate only one infrastructure without taking into account explicitly the interactions and interdependencies with other infrastructures. This paper proposes to develop a SINARI page 2/6 composite vulnerability model for the ICT or Cyber Infrastructure and the Power Systems Infrastructure using Complex Network Theory.

Complex Network is a new tool that allows modeling systems as graphs and has been extensively used to model, analyze, and understand large systems with non-trivial topologies and hidden interdependences [START_REF] Boccaletti | Complex Networks: Structure and Dynamics[END_REF]- [START_REF] Albert | Statistical mechanics of Complex Networks[END_REF]. This approach allows systems topology characteristics and connectivity properties to be known, as well as, fault and cascade phenomena analysis to be performed. This theory has been applied to analyze the vulnerability of Power Systems and Computer Networks [START_REF] Kröger | Vulnerable Systems[END_REF]. Critical Infrastructures can be modeled as a Complex Network of systems functioning together to achieve common purposes. In order to elaborate a joint model, this paper proposes to adapt the theory of Complex Numbers to the theory of Complex Networks. The result of this symbiosis is a twodimensional model, which allows inherent vulnerabilities of coupled infrastructures to be identified. Specifically, the proposed method bridges the Cyber Infrastructure and the Power Systems Infrastructure.

Section 2 presents an overview of Dependency and Vulnerability in Power Systems. Section 3 gives an introduction to Complex Networks and describes how the Complex Networks theory can be modified in order to introduce the Complex Numbers. Section 4 shows how this framework is applied to Power Systems. Finally, the conclusions are addressed in Section 5.

Dependency and Vulnerability

In order to identify, understand, and study vulnerabilities of Power Systems, it is important to know i) how the Power grid and the ICT inter-work; ii) the types of failures existing in both infrastructures; iii) the types of interdependencies; and iv) how these interdependencies affect the vulnerability of each infrastructure.

The main emphasis of Power Grids has been on providing a reliable, secure and economic supply of electrical energy to their customers. Control, protection, and monitoring systems are necessary for the Power Grid functioning. These systems are composed by ICT, which are defined as the technologies that involve acquiring, storing, processing, and distributing information by electronic means, including RTUs, IEDs, Computers, Servers, SCADA Systems, Routers, Gateways [START_REF] Bjorn | Overview of ICT components and its applications in electric power systems[END_REF]- [START_REF] Ericsson | Information security for Electric Power Utilities -results of Cigré WG D2.22[END_REF]. Therefore, there is a reciprocal relationship between both Infrastructures including two types of flows: energy and information. In one direction, Power Grid provides energy to ICTs; in the other direction, ICTs supervise, control, and manage the grid through commands and signals.

Three kinds of failures or outages can be found in Critical Infrastructures [START_REF] Rinaldi | Indentifying, understanding, and analyzing critical infrastructure interdependencies[END_REF]:

1. Cascading failures: Occur when a failure in one infrastructure causes a failure in a second infrastructure.

2. Escalating failures: Occur when a failure, resulting from the interaction between two infrastructures, exacerbates another failure.

3. Common-cause failures: Occur when two or more infrastructures are affected simultaneously because of an external and common cause, e.g. tornado, earthquake. Consequently, these failures show that the infrastructures are subjected to an increased risk from direct connectivity or spacial proximity, and that neighbor infrastructures are likely to be damaged after a single failure. In conclusion, Critical Infrastructures have strong interdependences, of which, there are four types: physical, cyber, geographic, and logical [START_REF] Rinaldi | Indentifying, understanding, and analyzing critical infrastructure interdependencies[END_REF].

Finally, how can all these dependencies and failures affect the vulnerability of an infrastructure? That is the main and still unanswered question. Coupled infrastructures vulnerabilities are the exploitable weaknesses facing a defined threat. Such threat, from a cyber-security point of view, can be either internal, i.e. people working in the organization and that have physical access rights; or external, e.g. DoS attacks.

Therefore, Power System and ICT Vulnerabilities have common threats that exploit the weakest element in the Coupled infrastructures. The weakest element is identified by its importance in the system. A node is "important" regarding its role in the network, as in [START_REF] Kröger | Critical infrastructures at risk: a need for a new conceptual approach and extended analytical tools[END_REF], the term 'importance' is intended to qualify the role that the presence and location of the node plays. There exist many methodologies in order to identify it, but they are made mostly for single infrastructures. However, our work with Complex Networks allows integrating both infrastructures in a single model, focusing on interdependencies and vulnerabilities that emerge from the interconnection of infrastructures.

Complex Networks

Methodology

A Network is a set of items with connections between them and a Graph is the network mathematical representation [START_REF] Newman | The Structure and function of complex networks[END_REF]. Networks are present in almost every aspect of life, e.g. traveling, calling by cell-phone, finding a job, chatting, etc. For this reason, mathematicians and experts in many domains have tried to find the best way to model real systems, considering the relations and dependencies between their components. Leonhard Euler solved the enigmatic seven bridges of Königsberg problem during the 18th century using what has been called the foundations of the Graph Theory.

Since its origin, the Graph Theory has been evolving thanks to the computerization of data acquisition and the availability of high computing power. Some of the main contributions include the Random Graphs Theory by Paul Erdös and Alfred Rényi in 1959 [START_REF] Erdös | On random graphs[END_REF], the introduction of the `Small-World' Concept by Watts and Strogatz (WS) [START_REF] Watts | Collective dynamics of 'Small-World' networks[END_REF], and the discovery of `scale-free' characteristic in large

i i i i A elec                              0 1 0 1 0 0 1 0 1 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 i i i i i i i i A ICT    V j
hj h a k networks by Albert Barabasi [16]. The application of Graph theory to study large and complex systems is called Complex Networks theory. A complete overview of Complex Networks history is presented in [START_REF] Watts | Collective dynamics of 'Small-World' networks[END_REF] and [START_REF] Wang | Complex networks: Smallworld, scale-free and beyond[END_REF].

FIG. 2 : Example System FIG. 3 : Example System Graph

Vertices and Edges compose Complex Networks. Vertices represent system elements such as buses, routers, airports or people. Edges represent the connections or relations between vertices; these connections can be physical, logical or functional. Some common edges include: power lines, optical fiber, flight itineraries, and friendship.

A graph G is a pair of sets (V,E). Where the elements of V ≡{v 1 , v 2 , v 3 ,…, v n } are the vertices or nodes and n is the number of vertices. While the elements of E ≡{e 1 , e 2 , e 3 ,…, e m } are the edges between the vertices and m is the number of edges.

A graph can be represented by a n×n matrix A, called adjacency matrix; where every row and column represents a vertex in the graph. Normally, its entry a hj is 1 if the edge exists between the hth and jth vertices and 0 otherwise. However, in order to create a Complex Network for two infrastructures, this paper proposes to set the entry a hj as a Complex Number (i.e. 1+1i). The matrix A is constructed according to (1), with w the number of outbound links from node h to node j, and x the number of inbound links from node h to node j [START_REF] Hoser | Eigenspectral Analysis of Hermitian Adjacency Matrices for the Analysis of Group Substructures[END_REF].

a hj = w + i x (1)
Since there are different layers in systems (electrical and ICT connections, logical and geographical dependencies), different adjacency matrices will represent each system. For example, Fig. 2 shows a basic electric grid and its ICT system. It is clear that there are many interdependencies between its components. For instance, Generation node supplies energy to loads, switches help to send and receive information from and to the control center, so electrical nodes depend on Control Centrer throught the ICT network. A single representation, as shown in Fig. 3, is not suitable to model all types of interdependencies. Therefore, this paper proposes to use multiple adjacency matrices. To show how to create these matrices, ( 2) and (3) show the adjacency matrices for the Electrical Layer and the ICT Layer shown in Fig. 4.

(

) (3) 2
Vulnerability assessment is a systematic evaluation in which quantitative or qualitative techniques are used to identify exploitable weaknesses in a system exposed to hazard. In order to do so, Power Grid and ICT infrastructures are modeled as Complex Networks using Complex Numbers. The assessed indexes are Node Degree and Efficiency. The higher the degree is (or lower the efficiency becomes), the more important the vertex in the system is. Then, an important vertex is highly vulnerable to a coordinated attack and/or a random failure.

FIG. 4 : Multi-Layer Analysis

The importance or prestige (for Social Networks) of a node is characterized by the node degree (number of inbounds and outbounds connections). According to [START_REF] Boccaletti | Complex Networks: Structure and Dynamics[END_REF], the degree k h of a node i is defined in terms of the adjacency matrix A as shown in (4). ( 4)

SINARI page 4/6      j h V j V h hj c e c e c d m m E , , 1 1 
    j h V j V h hj e c e c e d m m E , , 1 1
In the proposed method, the degree k h is a complex number (k h = y + iz) where y is the out-degree and z the indegree. Results from the test case are shown in Table 1.

TAB. 1 : Node Degree results

Vertex

v 1 v 2 v 3 v 4 v 5 v 6 k h Elec 1 2+i 1+i 0 i i k h ICT 0 1+i 1+i 1+i 3+3i 2+2i
Another method to identify important nodes is the Efficiency. The concept of Efficiency was introduced in [START_REF] Latora | Efficient behavior of smallworld networks[END_REF]. It is used to evaluate and measure how efficiently a node exchanges information with other nodes. The mathematical representation is in ( 5) and ( 6) for multiinfrastructures, where d hj is the shortest path length between h and j.

(5) ( 6) m e is the number of electrical edges, and m c is the number of ICT edges. For the test case, presented in Fig. 2, Table 2 presents the efficiency of nodes.

TAB. 2 : Efficiency results

Vertex

v 1 v 2 v 3 v 4 v 5 v 6 E e 0,

Test System and Results

Test System Description

The test system is a modification of the typical French Distribution Network presented in [START_REF] Stahl | Experimenting with infrastructures[END_REF]. The Power Grid has 14 power-bus, 17 lines, 7 distributed generation, 9 loads, and 3 transformers HTB / HTA, as shown in Fig. 5. Aside from this network, there is a considerable supporting ICT infrastructure which involves 2 routers from a communication network of a public ICT provider (nodes 19 and 20) and 1 WiMax BS (node 22). It includes as well multiple links ADSL, PSTN / ISDN, Optic Fiber, and Ethernet technologies. Router 23 represents the private LAN-GigaEthernet connecting the electrical buses 2, 3 and 4.

Complex Network Modeling

To create the complex network, the power-bus, loads, routers and Wimax are considered as vertices; and the power-lines, communication links (ADSL, Ethernet, etc.) are considered as edges. Fig. 6 shows the resulting graph. Electrical nodes are numbered from 1 to 14 and ICT nodes from 15 to 22.

FIG. 6 : Graph

Tables 3 and4 show the node degree and efficiency indices for the resulting graph. According to the tables nodes 2, 8 and 1 are the most critical for the Power System, and the failure consequences can be measured thanks to the node degree. For instance, to loose the node 2 means to loose 3 other nodes. However, it is more interesting the case of the set of nodes {8, 9, 19}, the three of them have a low efficiency index and the node 19 is supplied of energy from node 8, that depends on node 9 (source). Therefore, this method can help to identify important and critical elements on the coupled system, and to detect possible consequences.

TAB. 3 : Node Degree results

Vertex

k h Elec k h ICT 1 1 + 2i 0 2 3 + 2i 2 + 2i 3 1 + 1i 2 + 2i 4 1 + 1i 2 + 2i 5 1 + 1i 1 + 1i 6 1 + 1i 1 + 1i 7 1 1 + 1i 8 2 + 1i 1 + 1i 9 1 + 1i 1 + 1i 10 2

Conclusions

Vulnerability assessment for Critical infrastructures is an important task in order to guarantee the availability, reliability, and security. In order to do so, a new method to model Critical Infrastructures using Complex Networks have been presented. It shows a great flexibility and capability to perform vulnerability assessment for Power Systems.

The proposed method should serve to understand the interactions between different system components in a Power Systems, and to develop a risk analysis to identify ways to reduce the vulnerabilities. It will also help develop resilience studies, specifically in robustness, redundancy, and resourcefulness studies.

Modifications of this method can be made in order to model dependencies of different infrastructures, others than Power Grids and Cyber-systems. For instance, it can be applied to other Complex Networks parameters such as Centrality indexes, path length, clustering coefficient, or even to create dynamic multi-dimensional Complex Networks to model Cyber-attacks and their consequences. It could thus be extended to an n dimensional Complex Network to assess interactions between more than two infrastructures.

The significance of this paper lies in its ability to bridge two infrastructures in a single flexible model, which is an important step to study Critical Infrastructures vulnerability and to design safe and secure Power Systems, in this new era where utilities increasingly rely on the public Internet.

17 FIG. 5 :

 175 FIG. 5 : Power System and ICT networks


Axel Von Arnim

Mathias Perrollaz

Arnaud Bertrand

Jacques Ehrlich

Vehicle Identification Using Infrared Vision and Applications to Cooperative Perception

Vehicles will be in the next future equipped with V2V telecommunication means to exchange data, such as the presence of an obstacle on the road, or an emergency braking notification. Vehicles are also more and more equipped with perception systems (cameras, laser scanners, radars) that enable them to explore the immediate environment, including other vehicles. We propose in this paper an on-board optical vehicle identification system to enable telecom and perception systems to cooperate. The optical identification identifies which vehicle, in the scene captured by the perception system, is sending information via telecom.

I. INTRODUCTION

Using vehicle to vehicle telecoms to extend the perception of local on-board systems like vision [START_REF] Labayrade | Real time obstacle detection on non flat road geometry through 'v-disparity' representation[END_REF], laser scanners [START_REF] Kaempchen | Feature-level fusion for free-form object tracking using laserscanner and video[END_REF] or radars [START_REF] Takahashi | Obstacle detection using millimeter-wave radar and its visualization on image sequence[END_REF] is an attractive idea. For example (see Fig. 1), a collision mitigation application using stereovision techniques can take advantage of information about the braking power of the followed vehicle or an emergency brake event notification via telecoms [START_REF] Mourllion | A study on the safetycapacity tradeoff improvement by warning communications[END_REF]. But to be relevant in a complex road scene, where several vehicles are involved, the information has to be localized, otherwise we have the knowledge of a vehicle's behaviour (speed, braking action) but are unable to determine which vehicle of the scene is doing so. If the followed and the following vehicles are accurately and reliably localized, it is possible for the following vehicle to match, in the scene captured by perception sensors, which vehicle is sending relevant data. Unfortunately, GPS based localization is not accurate enough to discriminate between vehicles on different driving lanes, and if it were (differential GPS), it would not be reliable since the GPS signal is masked by tunnels, dense trees or urban canyons. So the need for an optical identification system, that delivers the localized identification of vehicles directly in the perception sensors captured scene is strong.

For now, the main systems developed use automatic license plate recognition [START_REF] Albiol | Automatic license plate reading using mathematical morphology[END_REF], or spatial pattern based techniques [START_REF] Bruce | Fast and accurate vision-based pattern detection and identification[END_REF]. These systems are limited by the accuracy of vision sensors (which induces a range limitation). Other systems use infrared light to extract the ID information from the view. But they use simple, often binary ID symbols to classify types of vehicles, rather than to individually identify them [START_REF] Usami | Development of cooperative vehicle following system for ahs in mixed traffic[END_REF]. A more suitable technology has been tested, that combines two IR beacons with a common signaling protocol by means of radio-communication to identify vehicles [START_REF] Takada | Multiple vehicle identification in a longitudinal ranging system[END_REF]. But this system is impeded by both the facts that it is Axel von Arnim (vonarnim@lcpc.fr), Mathias Perrollaz (mathias.perrollaz@inrets.fr) and Jacques Ehrlich are with LIVIC (INRETS / LCPC), Bât. 824, 14, route de la Minière 78000 Versailles -France. vulnerable to partial occultation, and that it needs a radio support signal to work. A system using two IR beacons at the rear of vehicles to measure inter-vehicle distances uses an optical identifier protocol (spatial and frequential), but the purpose is to match the two corresponding beacons on each vehicle [START_REF] Blmal | for Rockwell International Corporation. Inter-vehicle distance measuring system[END_REF]. So, a robust, vision-only, and long-range system does not, to our knowledge, exist. In this paper we present an optical identification sensor that enables middle to long-range identification (> 60m), does not rely on high resolution cameras and is not vulnerable to partial occultation. In part II we give a description of our sensor, from hardware to detailed software. Part III gives some experimental material. Finally, in part IV we present applications of the sensor to three cooperative perception systems.

II. OPTICAL VEHICLE IDENTIFICATION

A. General Overview of the Sensor 1) Working Principles: The goal is to identify vehicles, or more generally objects, in the road scene and localize them in the referential of the receiver-equipped vehicle. The need to identify vehicles in a camera captured scene, without having to detect spatial patterns, leads to a near IR solution. The proposed system is a cooperative sensor, composed of two parts: the emitter and the receiver.

Let us consider a car following another car and wanting to identify it. Both cars have telecom means.

On the leading car, we attach a led-based infrared emitter, which codes, using time-frequency digital coding (ie blinking) the unique identification number of the car, for example, its IP address. Time-frequency signals are more robustly and farther perceptible by a vision sensor than spatially coding signals, who are vulnerable to both partial occultation and vision sensors resolution. On the following vehicle, we use a low-resolution, high-frequency, infrared band pass filtered digital camera, which detects and decodes the signal, and localizes it in the same referential as the perception systems. Thus, each detected vehicle can be identified and the information from radio-communication and optical sources can be merged. An overview of the complete process, from image acquisition to delivery of the identified objects to other perception applications is presented on Fig. 2.

2) The emitter: The emitter (Fig. 3) is composed of a panel (named the lamp) of infrared emitting leds (SFH 484-2 emitting at peak 880 nm in our prototype), an electronic command controller that generates the coding signal, and a power controller delivering the necessary electrical power. Powerful, directional infrared LEDs can enable a longitudinal perception range of more than 300m, in good conditions. The digital coding signal conforms to a protocol (Fig. 4). A cycling frame is emitted permanently, containing a start sequence (4 bits), the usable code (identifier, as many bits as necessary), a 0 and a parity bit. 3) The receiver: The receiver is a high-frequency digital camera which captures the image of the road scene, and a processing algorithm which decodes and localizes the signal from the emitter. The camera is equipped with a bandpass infrared filter centered around 880 nm.

B. Low-level Detection Algorithm

A low-level algorithm (Fig. 2-a) processes the input image to extract the spots corresponding to emitters. Most false detections (other infrared sources like sky, vegetation) are eliminated by shape analysis. The order of operations is:

1) IR-filtered image is acquired. Emitters appear like bright round spots (Fig. 5-b). 2) Image is binarized using a threshold, so that only the brightest pixels are kept (the most IR emitting sources) (Fig. 5-c). 3) Regions are made more consistant, less fragmented, using morphological dilatation. 4) White pixels are clustered: each connected region is a labelled cluster (Fig. 5-d). 5) Regions are accepted or rejected depending on their size and roundness. Large regions, corresponding to the sky, and regions which don't look like a spot, for example pieces of vegetation, are rejected. The algorithm used to determine the roundness of a region computes the rate of pixels of the region inside its circumcircle. Accepted spots are named targets. Some objects, like solar reflections or vehicle or traffic lights appear like spots and are then accepted. We will see in the next paragraph that they are eliminated by the high-level tracks filtering algorithm.

C. High-Level Filtering, Tracking and Decoding

A high-level algorithm (Fig. 2-b) tracks detected targets over time (in successive images) to follow their movement and analyse their logical state (0=extinct, 1=lit). In the following, a tracked target is named a track.

The tracking is based on a very simple neighborhood prediction. The neighborhood in which a track is expected in the next image is a region around the last position where it was seen. As long as the track is extinct, the neighborhood is growing, faster horizontally than vertically, because the horizontal movement of the followed vehicle in the image is obviously greater than its vertical movement. When the track appears again, the neighborhood is reduced to the center of this new observation.

The observed states of each track are translated into bits, which are recorded into a FIFO. This resulting bitstream is continuously decoded as a message. Besides, a confidence is attached to each track. A new track starts with maximum confidence, and during the tracking, this confidence is managed as follows (Fig. 2-c):

• If the message does not conform to the communication protocol (Fig. 2-e), or if the parity bit of the message is incorrect (Fig. 2-f), the track either is a false detection or has a temporary communication error and is marked as incorrect. It will continue to be tracked. Its confidence keeps decreasing as long as the problem remains.

• If the track has been too long incorrect, ie. the confidence has reached zero, it is deleted and ignored in the next images.

• If the message conforms to the communication protocol and the parity bit is ok (Fig. 2-g), its code (identifier) and localization in the image are attached to it. If the track was incorrect before, its confidence increases, but it remains incorrect as long as the confidence has not reached its maximum value. The confidence increase is quicker than the decrease, to avoid accumulating too much delay when a track is temporarily incorrect due to short masking for example. Of course, The track will continue to be tracked in successive images. False tracks, created by solar reflections or vehicle and traffic lights that passed the low-level image analysis, are here eliminated because they don't blink, or they blink without conforming to the communication protocol (random blinking).

Finally, valid tracks are output (Fig. 2-h) to the absolute localization algorithm (Fig. 2-d) with their position in the image coordinate system and their identification number.

Figure 6 shows three tracks. One, in thick white, is a valid track and two, in thin grey, are invalid tracks. They correspond to the front lights of a vehicle coming on the opposite lane and are rightly rejected.

D. Absolute Localization

Once the emitter-equipped objects have been detected, identified and localized in the identification image, they have to be localized in the coordinate systems of the other sensors of the vehicle, in order to be shared (Fig. 2-d). For that purpose, an absolute coordinate system R a is defined. As all sensors are rigidly linked to the car frame, they can all refer to R a . Therefore, whichever sensor is used, the absolute localization problem finally consists in an inverse mapping of the identified tracks into R a . The position of the tracks can then be transferred into any sensor's coordinate system.

The only 2 camera coordinates (U id , V id ) are not sufficient to solve this 3 variable problem. This issue adds a new constraint on the emitter: its height Y L is supposed to be known. That constraint is not so strong, since the height of the emitter can be coded into the identifying code itself.

1) Geometrical Description:

The geometry of our receiver is presented on figure 7. The sensor is a video camera represented by a pin-hole model whose behaviors are: α u id and α v id , respectively the focal length measured in pixels horizontally and vertically (pixels are rectangular), and (U 0 id , V 0 id ) the projection of the optical center among the image plane. The camera is placed in R a on the (X 0 id , h id , Z 0 id ) position. Yaw and roll angles are considered null and constant. The pitch angle is θ id .

X a Z a -Y a h s h id X id Z id Y id θ id Z o id X o id Identification image X a Z a -Y a h s h id X id Z id Y id θ id Z o id X o id

Identification image

Fig. 7. Geometry of the identification sensor in the absolute coordinate system.

2) Coordinate Transform: The easiest approach to localization relies on a plane world assumption. Indeed, considering that the road surface is the (X a , 0, Z a ) plane, and according to the projection equations of the pin-hole model, the coordinates (X a , Y a =Y l , Z a ) of the emitter in R a can be easily computed using equation system 1.

       Z a = Z 0 id + (Y L +h id )(α v id cos θ id -(V id -V 0 id ) sin θ id ) (V id -V 0 id ) cos θ id +α v id sin θ id X a = X 0 id + (U id -U 0 id )((Y L +h id ) sin θ id +(Za-Z 0 id ) cos θ id ) α u id (1)
In operation, h id and θ id parameters can change significantly, due to vehicle pitch. Thus, a pitch sensor is required for a correct coordinate transform. For this purpose, any inertial sensor could be efficient. Stereovision is also well designed for that application. Indeed, as described in [START_REF] Labayrade | Real time obstacle detection on non flat road geometry through 'v-disparity' representation[END_REF], the height and pitch of the stereoscopic sensor can be retrieved using the v-disparity approach. h id and θ id are then immediately deduced, because the stereoscopic sensor is rigidly linked to the identification sensor.

The use of stereovision for the estimation of these parameters is an interesting solution because it could also be used in parallel for obstacle detection.

III. EXPERIMENTAL RESULTS

A. Scenarii

Here are the parameters which have an effect on the experiments:

• speed of receiver-equipped vehicle • camera iris aperture of the receiver. In all the following test scenarii, we use an auto-iris, to have optimal iris aperture for our application

• meteorological (light) conditions. A cloudy day has best light conditions, and a sunny day, the worst ones

• camera framerate of the receiver (related to exposure time), expressed in fps (frames per second)

• camera resolution of the receiver • emitter frequency, expressed in Hz, equivalent to bits/s. The ratio between the camera and the emitter frequencies gives the number of acquired images per bit emitted. This ratio, because of Shannon's theorem, must be greater than 2 Here are the features of the system that we wanted to evaluate. They are:

• static identification range: maximum range where a correct identification of the object is possible. Measured on static emitter-equipped object, with static receiver.

We placed an emitter-equipped vehicle at one side of a long straight portion of road, and we placed a receiver at different longitudinal positions on the road until we found the limit detection range.

• identification time: time between first detection of the emitter and correct identification. Measured on static emitter-equipped object, with static receiver. We chose an optimal position for the receiver.

• false detection rate: number of false detections per minute or km. Moving receiver (vehicle). We drove a receiver-equipped vehicle on public roads with normal to heavy traffic, during a reasonable time.

• non detection rate: number of emitter-equipped objects not detected. Moving receiver. Measured on the hereafter explained dynamic scenario.

• dynamic identification range: same feature as the static one, but measured on the dynamic scenario hereafter explained. Static emitter (road sign), moving receiver.

• localization accuracy: measured by placing an emitter on a vehicle moving away on a straight road. The results of the absolute localization (among Z a axis) are compared to a reference, measured using a lidar, with centimetric precision. Dynamic scenario on test track: we placed emitterequipped road signs on the roadside, driving a receiverequipped vehicle at various known speeds. Despite the system is mainly a vehicle identification sensor, we designed the scenario with static emitters, in order to control the previously listed parameters. And how could static emitters be better employed than being attached to road signs? But we also tested the system while following a car.

B. Results

In all experiments, we used a Pulnix greyscale high speed CCD camera, in binning mode (320x120) to increase the acquisition framerate and lower the CPU consumption. This low resolution is sufficient because we do not detect spatial patterns, but bright spots. The computer ran a Pentium D 940 with 1Gb RAM, on Windows XP. The emitters used 16 bit frames, including start code and parity bit.

1) Static experiments:

The main parameters here are a null speed of the receiver and the emitter (static), a variable camera framerate, thus a variable emitter frequency (according to Shannon's theorem). We tested the features with two values of the meteorological conditions parameters: cloudy and sunny day. The images per bit value is here a measured value (during a calibration procedure), and not the theoretical ratio between camera and emitter frequencies. The difference is due to discretization artefacts and threshold effects. We can see on table I that the receiver has an optimal identification range with an exposure time of 500µs, which gives the best contrast between emitters and environmental IR sources.

The theoretical identification time is given by the formula: t th id = n * ib/f c , where n is the total number of bits of a message, including parity bit and start code sequence, ib is the number of images per bit acquired, and f c is the camera framerate. The experimental identification time is always slightly greater than the theoretical one because the emitter and the receiver are not synchronized.

The measured identification time for our optimal configuration is 100ms, which is a value that enables to identify objects moving very fast, or to identify static objects while driving very fast by.

So for the next presented experiments, we used following combination of parameters: cameraf ramerate = 514f ps, emitterf requency = 210Hz. The receiver algorithm always runs in real-time, consuming an average 20% CPU.

2) Dynamic experiments:

Table II shows the false detection rate for a scenario on public road, by sunny weather (bad condition). We had a few false detections due to solar reflections, with very little duration. For other dynamic features, we used the formerly presented scenario with emitter-equipped road signs on our test track. The main parameters were: variable speed of receiver (from 50 to 110 km/h), cloudy day (good condition), and emitters' positions as shown on figure 8.

On table III, you can see the list of identified road signs: #12, 6, 1, 9, 7 and the average identification range when we drove along at different speeds. The standard deviation from the average value is very low, so we only report here the average values. This means that the tested features have similar values at all different speeds. The ranges depend on the location of the road signs. Road signs #12 and 7 are seen from far away (Fig. 8) because they are located on straight road portions. #1 appears on the image just after the tight turn, so the range is short.

Figure 9-top shows what we call a chronogram of the scenario and represents the identification code detected versus time. So you can see the history of the scenario: the time a road sign is detected, the duration of the correct detection (ie its identification range, knowing the speed of the vehicle). It shows how the id ranges and times are related to the positions of the road signs on the track. You can see that the vehicle is continuously detected and identified, apart from a short period of time in the tight turn, near T ime = 12. A short break in the identification of road sign #12 appears after T ime = 8, because the road sign is momentarily masked by the leading vehicle.

Other tests with this scenario were performed with less optimal meteorological conditions (sunny weather). The id ranges for road signs located on straight lines were reduced, but remained compatible with the scenario and let us drive the receiver at the same speeds, up to 110 km/h. Finally, we have similar results at different speeds of the receiver-equipped vehicle, and the ranges are relatively high, allowing interesting applications, as explained in part IV.

3) Localization Accuracy:

The performances of the localization algorithm have been assessed using stereovision as a height/pitch sensor. Figure 10 gives the results of the process among the Z a axis.

The localization algorithm gives satisfying results on short ranges (until 25 meters). For longer ranges, the accuracy decreases because of the low vertical accuracy of each pixel. As expected, the accuracy is better when the identification camera is configured in binning mode 2-2 (320*240 pixels) than in binning mode 2-4 (320*120 pixels), but with no spectacular improvement.

IV. APPLICATIONS TO COOPERATIVE PERCEPTION

A. Telecom-Enriched Cooperative Detection

Let us consider a vehicle equipped with an obstacle detection system, telecommunication means and our identification system to combine them. As presented in figure 11, the localized identifiers can be associated with the corresponding detected objects. Then, the telecom messages can be matched with their detected senders.

The vehicle gets low level info from the detection system (relative position and speed), while telecommunications can provide high level information such as event notifications, masked obstacle alerts or characteristics of vehicles. Thus, as we have multilevel info, the vehicle can have a more adapted and reactive response to a change of its environnement. This enriched perception system is adapted to design a collision anticipation application. By extension, when applying this system to a pool of vehicles in a driving lane, it is possible to dramatically reduce the global time of reaction on a road event. This kind of strategy can strongly reduce the gravity of accidents, as reviewed in [START_REF] Mourllion | A study on the safetycapacity tradeoff improvement by warning communications[END_REF].

B. Cooperative Road Sign

If an emitter is placed on a road sign, you can detect and identify it with the onboard receiver in vehicles. So the application is to automatically detect and identify road signs or any other roadside object in the vehicle. For example, a vehicle detecting and identifying a dangerous curve road sign can warn the driver, or automatically reduce its speed. This idea is possible with any kind of road sign.

C. Long Range Cooperative Detection

As our optical identification sensor is able to detect equipped vehicles up to 300 meters, it can be used for long range obstacles detection. A basic approach consists in just localizing the identified tracks in R a . This kind of detection is not precise enough for driving assistance purposes, but provides seeds for more accurate vision systems. For example, when stereovision is also in use in the vehicle, the identified tracks, localized in R a can be used to define a region of interest in the stereoscopic images. Then, the distance to the obstacle is more precisely measured by stereovision. Moreover, using a confirmation strategy [START_REF] Perrollaz | Long range obstacle detection using laser scanner and stereovision[END_REF] can lead to a very robust detection, extending or replacing lidar capabilities.

V. CONCLUSION

This paper has presented a vision-only identification system, which is long-range and computation cost effective. It has proved to be efficient in many road scenarii. It enables to combine V2V telecoms with shorter range perception systems to make possible applications like extended or enriched perception. Besides, other applications like smart road signs have been successfully tested. Future works will deal with a complete evaluation of the performances of the sensor and the development of combined applications with V2V telecoms.

Fig. 1 .

 1 Fig. 1. Combining perception and telecoms.

Fig. 2 .

 2 Fig. 2. Overview of the identification algorithm.

Fig. 3 .

 3 Fig. 3. Prototype of the emitter lamp.

  Code '0' parity . . .

Fig. 4 .

 4 Fig. 4. Emitter signal protocol.

Fig. 5 .

 5 Fig. 5. Low-level Image Analysis.

Fig. 6 .

 6 Fig. 6. Two bad tracks (thin) and one valid (thick).

Fig. 8 .

 8 Fig. 8. Experimental Scenario.

Figure 8

 8 Figure 8 shows the test track, and the positions of the emitter-equipped road signs. Each road sign was emitting an identifying code (number) which is shown on the figure. The various positions of the emitters ensured that the receiver detected them with various light expositions.

Fig. 9 .

 9 Fig. 9. Chronograms of the scenario at speed = 90 km/h (top), and following a vehicle (bottom).

Figure 9 -

 9 Figure 9-bottom is another chronogram of the same scenario, but this time we followed an emitter-equipped vehicle. It is intended to show the quality and continuity of the leading vehicle identification.

Fig. 10 .

 10 Fig. 10. Results of the absolute localization step in Za, compared to a reference measured with a lidar.

Fig. 11 .

 11 Fig. 11. Telecom enriched detection: a) a car is detected thanks to the lidar and a message is received from vehicle 8. b) an object identified as 8 is perceived and localized in the scene. c) both detections are associated: the car in front is sending an emergency brake message.

TABLE I STATIC

 I EVALUATION OF THE PERFORMANCES.

	Camera Frame rate	410	514	595	650
	Camera exposure time (µs)	1000	500	250	100
	Emitter frequency (Hz)	175	210	250	275
	Images per bit	2.29	2.36	2.32 2.17
	Id. range -Best case (m)	360	400	378	83
	Id. range -Worst case (m)	65	110	107	62
	Theoretical min Id. time (ms)	90	73	62	53
	Id. time (ms)	105	100	76	65

TABLE II FALSE

 II DETECTIONS ON PUBLIC ROADS (WORST CASE).

	Total scenario duration	16 min
	Total scenario length	21 km
	False detection occurrences	5
	Total false detections duration	138 ms
	Average false detection duration	28 ms
	False detection rate (nb/minute)	0.31
	False detection rate (nb/km)	0.24

VI. ACKNOWLEDGMENTS

The authors particularly want to thank J.M. Chevreau and A. Douence for their precious help in the development of the prototypes of the system.

This study was partly realized in the European TRACKSS project, and a patent (fr-07-52663) has been filled.


Muhammad Marwan

Muhammad Fuad

Pierre-François Marteau
email: pierre-francois.marteau@univ-ubs.fr

TOWARDS A FASTER SYMBOLIC AGGREGATE APPROXIMATION METHOD

Keywords: Time Series Information Retrieval, Symbolic Aggregate Approximation, Fast SAX

The similarity search problem is one of the main problems in time series data mining. Traditionally, this problem was tackled by sequentially comparing the given query against all the time series in the database, and returning all the time series that are within a predetermined threshold of that query. But the large size and the high dimensionality of time series databases that are in use nowadays make that scenario inefficient. There are many representation techniques that aim at reducing the dimensionality of time series so that the search can be handled faster at a lower-dimensional space level. The symbolic aggregate approximation (SAX) is one of the most competitive methods in the literature. In this paper we present a new method that improves the performance of SAX by adding to it another exclusion condition that increases the exclusion power. This method is based on using two representations of the time series: one of SAX and the other is based on an optimal approximation of the time series. Pre-computed distances are calculated and stored offline to be used online to exclude a wide range of the search space using two exclusion conditions. We conduct experiments which show that the new method is faster than SAX.

INTRODUCTION

Similarity search is a hot topic in computer science. This problem has many applications in multimedia databases, bioinformatics, pattern recognition, text mining, computer vision, and others. Small structured databases can be handled easily. But managing large unstructured, or weakly structured, databases requires serious effort, especially when they contain complex data types.

Time series are data types that appear in many applications. Time series data mining includes many tasks such as classification, clustering, similarity search, motif discovery, anomaly detection, and others. Research in time series data mining has focussed on two aspects; the first aspect is the dimensionality reduction techniques that can represent the time series efficiently and effectively at lower-dimensional spaces. Different indexing structures are used to handle time series. Time series are high dimensional data , so even indexing structures can fail in handling these data because of what is known as the "dimensionality curse" phenomenon. One of the best solutions to deal with this phenomenon is to utilize a dimensionality reduction method to reduce the dimensionality of the time series, then to use a suitable indexing structure on the reduced space.

There have been different suggestions to represent time series in lower dimensional spaces. To mention a few; Discrete Fourier Transform (DFT) [START_REF] Agrawal | Efficient Similarity Search in Sequence Databases[END_REF] and [START_REF] Agrawal | Fast Similarity Search in the Presence of Noise, Scaling, and Translation in Time-series Databases[END_REF], Discrete Wavelet Transform (DWT) [START_REF] Chan | Efficient Time Series Matching by Wavelets[END_REF], Singular Value Decomposition (SVD) [START_REF] Korn | Efficiently Supporting ad hoc Queries in Large Datasets of Time Sequences[END_REF], Adaptive Piecewise Constant Approximation (APCA) [START_REF] Keogh | Locally Adaptive Dimensionality Reduction for Similarity Search in Large Time Series Databases[END_REF], Piecewise Aggregate Approximation (PAA) [START_REF] Keogh | Dimensionality Reduction for Fast Similarity Search in Large Time Series Databases[END_REF] and [START_REF] Yi | Fast Time Sequence Indexing for Arbitrary Lp norms[END_REF], Piecewise Linear Approximation (PLA) [START_REF] Morinaka | The L-index: An Indexing Structure for Efficient Subsequence Matching in Time Sequence Databases[END_REF], Chebyshev Polynomials (CP) [START_REF] Cai | Indexing Spatio-temporal Trajectories with Chebyshev Polynomials[END_REF].

Among the different representation methods, symbolic representation has several advantages, because it allows researchers to benefit from textretrieval algorithms and techniques that are widely used in the text mining and bioinformatics communities [START_REF] Keogh | Locally Adaptive Dimensionality Reduction for Similarity Search in Large Time Series Databases[END_REF].

The other aspect of research in time series data mining is similarity distances. There are quite a large number of similarity distances; some of them are applied to a particular data type, while others can be applied to different data types.

Among the different similarity distances, there are those that can be used on symbolic data types. At first they were available for data types whose representation is naturally symbolic (DNA and proteins sequences, textual data…etc). But later these symbolic similarity distances were extended to apply to other data types that can be transformed into sequences by using an appropriate symbolic representation technique.

In time series data mining there are several symbolic representation methods. Of all these methods, the symbolic aggregate approximation method (SAX) [START_REF] Lin | A Symbolic Representation of Time Series, with Implications for Streaming Algorithms[END_REF]) stands out as one of the most powerful methods. The main feature of this method is that the similarity distance that it uses is easy to compute, because it uses statistical lookup tables.

In this paper we present a new method that speeds up SAX by adding a new exclusion condition that increases the exclusion power of SAX. Our method keeps the original features of SAX, mainly its speed, since the new exclusion condition is precomputed offline.

The rest of this paper is organized as follows: in section 2 we present background on dimensionality reduction, and on symbolic representation of time series in general, and SAX in particular. The proposed method is presented in section 3. In section 4 we present some of the results of the different experiments we conducted. The conclusion is presented in section 5.

BACKGROUND

Representation Methods

Managing high-dimensional data is a difficult problem in time series data mining. Representation methods try to overcome this problem by embedding the time series of the original space into a lower dimensional space. Time series are highly correlated data, so representation methods that aim at reducing dimensionality by projecting the original data onto lower dimensional spaces and processing the query in those reduced spaces is a scheme that is widely used in time series data mining community.

When embedding the original space into a lower dimensional space and performing the similarity query in the transformed space, two main sideeffects may be encountered; false alarms, also called false positivity, and false dismissals. False alarms are data objects that belong to the response set in the transformed space, but do not belong to the response set in the original space. False dismissals are data objects that the search algorithm excluded in the transformed space, although they are answers to the query in the original space. Generally, false alarms are more tolerated than false dismissals, because a post-processing scan is usually performed on the results of the query in the transformed space to filter out these data objects that are not valid answers to the query in the original space. However, false alarms can slow down the search time if the algorithm returns too many of them. False dismissals are a more serious problem and they need more sophisticated procedures to avoid them.

False alarms and false dismissals are dependent on the transformation used in the embedding. If f is a transformation from the original space ) , ( in order to guarantee no false dismissals this transform should satisfy:

) , ( )) ( ), ( ( 2 1 2 1 u u d u f u f d orig trans ≤ , orig S u u ∈ ∀ 2 1 , (1)
The above condition is known as the lower-bounding lemma. [START_REF] Yi | Fast Time Sequence Indexing for Arbitrary Lp norms[END_REF] If a transformation can make the two above distances equal for all the data objects in the original space, then similarity search produces no false alarms or false dismissals. Unfortunately, such an ideal transformation is very hard to find. Yet, we try to make the above distances as close as possible.

The above condition can be written as:

1 ) , ( )) ( ), ( ( 0 2 1 2 1 ≤ ≤ u u d u f u f d orig trans (2)
A tight transformation is one that makes the above ratio as close as possible to 1.

SAX

Symbolic representation of time series has attracted much attention recently, because by using this method we can not only reduce the dimensionality of time series, but also benefit from the numerous algorithms used in bioinformatics and text data mining. However, first symbolic representation methods were ad hoc and did not give satisfactory results. But later more sophisticated methods emerged. Of all these method, the symbolic aggregate approximation method (SAX) is one of the most powerful ones in time series data mining.

SAX is based on the fact that normalized time series have highly Gaussian distribution [START_REF] Larsen | An Introduction to Mathematical Statistics and Its Applications[END_REF], so by determining the breakpoints that correspond to the chosen alphabet size, one can obtain equal sized areas under the Gaussian curve.

SAX is applied in the following steps: in the first step all the time series in the database are normalized. In the second step, the dimensionality of the time series is reduced by using the PAA. In PAA the times series is divided into equal sized segments or "frames" and the mean value of the points that lie within the frame is computed. The lower dimensional vector of the original time series is the vector whose components are the means of all successive frames of the time series. In the third step, the PAA representation of the time series is discretized. This is achieved by determining the number and the location of the breakpoints. This number is related to the desired alphabet size (which is chosen by the user), i.e. alphabet_size=number(breakpoints)+1

. Their locations are determined by statistical lookup tables, so that these breakpoints produce equal-sized areas under the Gaussian curve. The interval between two successive breakpoints is assigned to a symbol of the alphabet, and each segment of the PAA that lies within that interval is discretized by that symbol.

The last step of SAX is using the following similarity distance:

2 1 )) , ( ( ) , ( ∑ = ≡ N i i i t s dist N n T S MINDIST ( 3
)
Where n is the length of the original time series, N is the number of the frames, S ~and T ~are the symbolic dist is implemented by using the appropriate lookup table.

We need to mention that the similarity distance used in PAA is:

∑ = - = N i y x N n Y X d 1 2 ) ( ) , ( (4)
Where n is the length of the time series, N is the number of frames. It is proven in [START_REF] Keogh | Dimensionality Reduction for Fast Similarity Search in Large Time Series Databases[END_REF] and [START_REF] Yi | Fast Time Sequence Indexing for Arbitrary Lp norms[END_REF] that the above similarity distance is a lower bounding of the Euclidean distance applied in the original space of time series . This results in the fact that MINDIST is also a lower bounding of the Euclidean distance, because it is a lower bounding of the similarity distance used in PAA. This guarantees no false dismissals. Figure .1 illustrates the different steps of SAX.

EXPERIMENTS

We conducted experiments on different datasets available at UCR (UCR Time Series datasets) and for all alphabet sizes, which vary between 3 (the least possible size that was used to test the original SAX) to 20 (the largest possible alphabet size). We compared the speed of our method FAST_SAX, with that of SAX as a standalone method. The comparison was based on the number of operations that each method uses to perform the similarity search query. Since different operations take different execution times, we used the concept of latency time [START_REF] Schulte | Performance Evaluation of Decimal Floating-point Arithmetic in IBM Austin Center for Advanced Conference[END_REF]. We report in Tables 1 and in Figure 2 the results of (wafer). We chose to present the results of this dataset because it is the largest dataset in the repository. Also it is shown in (Muhammad Fuad and Marteau 2008) that the best results obtained with SAX were with this dataset.

The codes of SAX were optimized versions of the original codes, since the original codes were not optimized for speed The results of testing our method on other datasets are similar to the results obtained with (wafer)

The results shown here are for alphabet size 3 (the smallest alphabet size), 10 (the largest alphabet size in the first version of SAX), and 20 (the largest alphabet size in the second version of SAX)

The results obtained show that FAST_SAX outperforms SAX for the different values of i and for the different values of the alphabet size.

CONCLUSION AND PERSPECTIVES

In this paper we presented a method that speeds up the symbolic aggregate approximation (SAX). We conducted several experiments of times series similarity search, with different values of the alphabet size and different threshold values. The results obtained show that the new method is faster than the original SAX The future work can focus on extending this method to other representation methods, other than SAX, also on using other approximating functions.

Figure 1 :

 1 Figure 1: The different steps of applying SAX

Figure 2 :

 2 Figure 2: Comparison of the latency time between SAX and FAST_SAX for alphabet size=3 (a), alphabet size=10 (b) and alphabet size=20 (c)

Table 1 :

 1 Comparison of the latency time between SAX

	and FAST_SAX for i=1:4 and alphabet size=3, 10, 20
		g= 3	g= 10	g= 20
	FAST_SAX 1.0592E6 3.7136E5 2.6734E5
	SAX	5.1291E6 1.5764E6 1.2759E6
		i=1
		(a)
		g= 3	g= 10	g= 20
	FAST_SAX 7.2062E6 3.3509E6 2.2255E6
	SAX	1.567E7	4.8078E6 3.4253E6
		i=2
		(b)
		g= 3	g= 10	g= 20
	FAST_SAX 1.3717E7 1.1444E7 9.5446E6
	SAX	2.1944E7 1.4428E7 1.1144E7
		i=3
		(c)
		g= 3	g= 10	g= 20
	FAST_SAX 2.2697E7 1.6928E7 1.6611E7
	SAX	2.8287E7 2.2179E7 1.9877E7
		i=4
		(d)

THE PROPSED METHOD

The basis of our method is to improve the performance of SAX by combining it with an exclusion condition that increases its pruning power, and by using different reduced spaces, and not only one.

We divide each time series into N segments. Each segment is approximated by a polynomial. In this paper we use a polynomial of the first-degree for its simplicity, but other approximating functions can be used as well.

Since this approximating function is the optimal approximation of the corresponding segment, the distance between this segment and this approximating function is minimal. The time series is also represented using SAX in the way indicated in section 2.2. A polynomial of the same degree is used to approximate all the segments of all the time series in the database. So now each time series has two representations: the first is a n -dimensional one, by using the approximating function, and the second is a N-dimensional one, by using SAX. We also have two similarity distances; the first one is the Euclidean distance, which is the distance between a time series and its approximating function, and the other one is MINDIST, which is the distance in the reduced space.

Given a query ) , ( ε q , let u , q be the projections of the u , q , respectively, on their approximating functions, where u is a time series in the database.

By applying the triangular inequality we get:

Taking into consideration that u is the best approximation of u , we get:

By substituting the above relation in ( 5), we can safely exclude all the times series that satisfy:

In a similar manner, and since q is the best approximation of q , we can safely exclude all the time series that satisfy:

Both ( 7), ( 8) can be expressed in one relation:

In addition to the exclusion condition in (9), since MINDIST is lower bounding of the original Euclidean distance, all the time series that satisfy:

Should also be excluded, Relation (10) defines the other exclusion condition.

The Offline Phase: The application of our method starts by choosing the lengths of segments. We associate each length with a level of representation.

The shortest lengths correspond to the lowest level, and the longest lengths with the highest levels. Each series in the database is represented by a first-degree polynomial, which is the same approximating function for all the time series in the databases. The distances between the time series and their approximating function are computed and stored. In order to represent the time series, we choose the alphabet size to be used. SAX appeared in two versions; in the first one the alphabet size varied in the interval (3:10), and in the second one the alphabet size varied in the interval (3:20). We choose the appropriate alphabet size for this datasets. The time series in the database are represented using SAX on every representation level

The Online Phase: The range query is represented using the same scheme that was used to represent the time series in the database. We start with the lowest level and try to exclude the first time series using (9) if this time series is excluded, we move to the next time series, if not, we try to exclude this time series using relation (10). If all the time series in the database have been excluded the algorithm terminates, if not, we move to a higher level. Finally, after all levels have been exploited, we get a potential answer set, which is linearly scanned to filter out all the false alarms and get the true answer set.


Blandine Romain

Veronique Letort

Olivier Lucidarme

Florence D'alché-Buc

Laurence Rouet

Véronique Letort

Florence D'alché-Buc

Registration of free-breathing abdominal 3D contrast-enhanced CT

Keywords: free-breathing, spatio-temporal registration, DCE-CT, evaluation criteria

  L'archive ouverte pluridisciplinaire

Introduction

Functional imaging has gained attention for oncology therapy. It is based on the acquisition of time sequences combined with contrast injection. Long acquisitions (up to 5 minutes) allow to study the complete dynamics of diffusion, including the effects of tissue permeability. However, it is possible to get relevant information only if data are correctly registered. Two main approaches are possible to deal with respiratory motion, sequential breath hold [START_REF] Koh | Dynamic contrast-enhanced ct imaging of hepatocellular carcinoma in cirrhosis: feasibility of a prolonged dual-phase imaging protocol with tracer kinetics modeling[END_REF] or free-breathing. In the present study, free breathing is selected, since a less stressed patient usually provides sequences with more regular motion. The current study focuses on tumors in liver or kidney using dynamic enhanced CT. Methods for registration of images acquired with free-breathing are being developed essentially for contrast-enhanced MRI. Wollny et al. [START_REF] Wollny | Exploiting quasiperiodicity in motion correction of free-breathing myocardial perfusion[END_REF] exploit the quasi periodicity within myocardial perfusion and Li et al. [START_REF] Li | 3d non-rigid motion correction of free-breathing abdominal dce-mri data[END_REF] explicitly estimate the respiratory cycle within abdominal perfusion. In both cases, an important assumption is a high time resolution.

In CT, it is also possible to have high time resolution, but this is combined with a short time of acquisition (2 minutes) [START_REF] Romain | Registration and functional analysis of ct dynamic image sequences for the follow-up of patients with hepatic tumors undergoing antiangiogenic therapy[END_REF], and this does not allow to extract permeability. For CT perfusion, with the objective of extracting adequate data for functional analysis, a longer time of acquisition is required which implies, mostly for dose considerations, a low time resolution. The second constraint is a limited imaging field of view in z-direction. To the best of our knowledge, no previous studies have investigated image registration on free-breathing with abdominal dynamic contrast-enhanced CT (DCE-CT) for long acquisitions. In that context, the aim of this paper is, in a first part, to propose the complete pipeline of the intensity-based methods that we have proposed to compensate free respiratory motion. The second main contribution of this study is the definition and evaluation of a dedicated similarity measure adapted from the difference entropy (DE), which is compared to the classically used DE [START_REF] Rodrigues | Evaluation of similarity measures in contrast enhanced echocardiography motion detection and registration[END_REF]. Finally, we propose the definitions of a temporal and a spatial criteria, in order to objectively evaluate the quality of the registration and to compare the results obtained when testing the different similarity measures.

MATERIALS AND METHODS

Data acquisition

39 data acquisitions were performed with a 256 slices CT (Brilliance iCT 256, Philips Healthcare, The Netherlands) on 14 patients. All patients gave informed consent. The dynamic CT protocol consists in acquiring 48 volumes every 2.5 seconds then 18 volumes every 10 seconds (80 kV, 80 mAS, rotation time of 0.33 seconds, dose of iodine 80 cc). With 8 cm detector coverage in z-direction, the CT scanner allows for 5.5 cm effective z-coverage in a single rotation with a 3D axial cone beam correction. This coverage is large enough to keep most of the lesion in image volume. In the axial plane, data are reconstructed with a pixel size of 0.68 mm × 0.68 mm [START_REF] Romain | Optimisation of ct reconstruction for the registration of ct liver perfusion sequences[END_REF]. Reconstructed volumes are expressed in Hounsfield Units (HU). Image intensities can thus range from -1024 (air) to 2000 HU (bone).

Spatio-Temporal registration

The main challenges for registration in these conditions are:

-Limited volume height: data acquisition covers about 6 cm in z-direction, while the amplitude of respiratory motion may reach 2 to 3 cm, in the same z-direction. This means that the whole liver or kidneys are not imaged. Consequently, full organ tracking is not possible and local approaches have to be used. -Contrast agent injection: intensities of a given tissue vary with time (figure 1). The choice of a similarity measure which is robust to contrast-induced intensity variations is essential in our method.

-Free breathing and low time resolution: since patients breath freely, a first consequence is that acquisitions correspond to different instants in respiratory cycle and a prediction using a respiratory model is not possible. A second consequence is the existence of strong intensity variations between two consecutive volumes due to contrast arrival and low time resolution (2.5 or 10 seconds). In such a context, the problem of intensity differences is present even between successive pairs of volumes (I t onto I t-1 ) and registering these pairs of volumes would lead to unnecessary errors accumulation. So, it is preferable to align all volumes onto a common reference volume I r . Since the majority of image volumes are enhanced, an intuitive solution would be to choose a reference volume I r with contrast (t ≈ 30). However, the assumption that the tumor is well centered can only be done in the first volume, which is why we selected it as the reference volume I r . In order to deal with these challenges, we present an adapted spatio-temporal registration pipeline. Each reconstruction provides a dynamic sequence DS I = (I t , t ∈ (1, ...T )), where I t is a 3D image acquired at time t and T is the number of acquisition times. Our registration approach involves two main steps, followed by a step of regularization and warping:

1. Global z-translation : since the main motion is in the z-direction, the first step consists in a global evaluation of the z-translation. To improve this initial step, thresholds are applied to exclude air (∼-1024 HU) and intense contrasts (≥600 HU), such as the contrast of the static aorta in z-direction. We considerer all pixels which intensity is between -800 HU and 600 HU.

Since we are considering one global translation for the whole masked volume at this point, we verified that the local contrast diffusion was not an important disturbance, and that the sum of square of differences (SSD) was the most appropriate metric (compared to entropy and mutual information, results not shown).

2. Multi-resolution blockmatching [START_REF] Tzovaras | Evaluation of multiresolution block matching techniques for motion and disparity estimation[END_REF]: initialized by the z-translation found in the first step, this 3D registration method is computed with block sizes of (11 × 11 × 7 mm 3 ). It provides, for each block, motion vectors corresponding to translations in x, y and z directions. 3. Regularization: a regularization step is necessary to smooth the motion vector fields. A gaussian filter is computed in 3D (σ = 1.2). 4. Warping: warping consists in reconstructing new image sequences (J 1 , ...J T )

with respiratory motion compensation. A trilinear interpolation is used.

Masking

In order to improve the algorithm robustness and its computing efficiency, masked volumes are computed. For each volume, a mask is defined by excluding spine, which have a different motion from the rest of the abdomen, and background, which does not give useful information.

Our mask consists of two 3D zones with z-translation invariance. These two zones can therefore be represented in an axial plane of arbitrary z-value.

On each axial plane, a threshold of 500 HU is first applied to extract the bones: the largest connected component is the spine. Similarly, the largest connected component of pixels of intensity values lower than -500 HU are considered to belong to the background (figure 2).

Fig. 2. Mask example: white zone is the valid region, excluding spine and background.

Similarity measures for multi-resolution blockmatching

To locally register the volumes I r and I t , our approach is based on a similarity measure which is potentially robust to intensity variations, namely DE. The entropy H of a discrete random variable X with probability mass function p(X) is defined as:

H(X) = E(-log p(X)) (1)
The difference between I r and I t , noted I r-t , is used to evaluate DE. The classical definition is based on the normalized histogram hist(I r-t ) of I r-t with B bins:

DE = - B j=1 hist(I r-t ) j * log(hist(I r-t ) j ) (2)
DE is classically used for registration, in particular in contrast enhanced echocardiography [START_REF] Rodrigues | Evaluation of similarity measures in contrast enhanced echocardiography motion detection and registration[END_REF].

However, it is known that the quality of registration is very dependent on the number of bins B. So, our approach is to approximate DE to have a more dedicated and robust similarity measure. We can define H directly in the image domain Ω as:

H = 1 card(Ω) card(Ω) i=1 -log p i ( 3
)
with p i the intensity of pixel i and card(Ω) the number of pixels Ω.

In addition, assuming that the pixel intensities in each block are distributed following an univariate Gaussian distribution with mean µ and variance σ, DE G is expressed as:

DE G = log(σ) + 1 card(Ω) card(Ω) i=1 (p i -µ) 2 2σ 2 (4)
We can then compare our proposed similarity measure DE G and the classical one DE. To this end, we define two evaluation criteria.

Evaluation criteria

In order to evaluate the quality of registration for each similarity measure, we define two criteria, related to the temporal and spatial dimensions of the method.

Temporal evaluation criterion:

We propose a new evaluation criterion based on curve smoothness. The quality of the registration is assessed by the smoothness of the time curve of mean intensities after registration (figure 3). Indeed, after registration, time-intensity curves should reflect the contrast intake only, and not the respiratory motion. Since the presence of kidney in a given axial slice is highly sensitive to respiratory motion, a region of interest (ROI) is defined within the kidney. Inside this ROI, mean intensities Īt of any image (I t ) before registration and Jt of the same image after registration are calculated, for each acquisition time t and plotted, see figure 3 for an illustrative example. The two curves are independently filtered by a median time-filter with a window span of 3 points, noted M (I t , 3) and M (J t , 3). Therefore, the sum of the absolute difference between the smoothed curve and the original curve can be used as indicator of the smoothness of the enhancement curve. Thus, we define SC as:

SC = T t=1 | Jt -M (J t , 3)| - T t=1 | Īt -M (I t , 3)| T t=1 | Īt -M (I t , 3)| (5)
The higher SC is, the better quality of registration is.

Spatial evaluation criterion: the usual metric of SSD was not used in our study since a change in contrast intensity would fool the evaluation. In order to remove variations due to contrast arrival, we decided to used the gradient of the volumes of differences. In order to calculate a spatial criterion on the quality of registrations, we computed two gradient image volumes G It-Ir and G Jt-Ir , and their associated norm ||G It-Ir || and ||G Jt-Ir ||. Eventually, we compute this spatial evaluation criterion, named Gradient Criterion (GC), as:

GC = 1 T T t=1 ||G It-Ir || -||G Jt-Ir || ||G It-Ir || (6)
G I is defined as the classical gradient of an image I:

G I = ∇ X (I) 2 + ∇ Y (I) 2
where ∇ X (I) is the gradient of I in x direction and ∇ Y (I) is the gradient of I in y direction.

From this gradient image, the L 2 norm is derived:

||G I || = (x,y,z)∈Ω G I (x, y, z) 2
This definition is based on the idea that, after registration, the values of difference volume gradients should be lower than before registration. So, a high GC shows a better quality of registration. Figure 4 shows the results of registration using a checkerboard. After registration, frontiers of organs are better aligned. After global z-translation, we apply our algorithm of multiresolution blockmatching with DE and our modified similarity measure DE G . In figure 5, results show that our proposed measure DE G is better adapted as similarity measure, according to SC and GC criteria. Note that we tested the assumption of gaussian distribution of pixel intensities in each block. It obviously depends on the size of blocks and their location (blocks that cover two different regions of the image have less normal distribution than blocks that include only homogeneous regions). Our algorithm was nevertheless robust even for blocks that did not verified this assumption.

CONCLUSION

We have presented a full setup to provide automatic registration of 4D DCE-CT sequences, using an algorithm that consists of global z-translation and multiresolution blockmatching. One main challenge, related to contrast arrival, has been specifically evaluated in the step of multi-resolution block-matching. In this part, we presented an adapted similarity measure (DEG), which was found to be more adapted to registration on free-breathing abdominal 4D contrast enhanced CT than DE.

Our mid-term objective is to provide a prototype that could be used by clinicians. Therefore, the computation time was one of our main constraints. Our method is relatively efficient. To register one image volume (512*512*22), it takes about 3 seconds (Intel(R) Core(TM) i5-2450M CPU, 2,5 GHz). This first result is encouraging for a future clinical use. A possible further refinement of our work could be to consider changing of similarity measures with time, according to the diffusion phases of the contrast.

Fig. 1 .

 1 Fig. 1. Image data at time 0 (left) and 15 seconds after injection (right). On the right figure: kidney is enhanced (white square) and aorta with hepatic artery are enhanced (white rectangle, middle).

Fig. 3 .

 3 Fig. 3. Mean intensities inside ROI on an axial slice during time of acquisition before registration Īt (dash line) and after registration Jt (solid line).

Fig. 4 .

 4 Fig. 4. On the top: before registration; on the bottom: after registration. At left: Checkerboards between image of reference (yellow) and current image (gray). At middle: zoom on kidney, which is better aligned after registration than before, especially in the red squares. At right, zoom on jejunum

Fig. 5 .

 5 Fig. 5. Values of the two criteria, namely intensity curve smoothness SC (left) and gradient criterion GC (right) for the two similarity measures DE and DEG. These plots show that the proposed similarity measure DEG improves the evaluation criteria on the all 39 data acquisitions.


Frédérique Billy

Jean Clairambault

Olivier Fercoq

Tommaso Lorenzi

Alexander Lorz

Benoît Perthame

Modelling targets for anticancer drug control optimisation in physiologically structured cell population models

Keywords: Cell population dynamics, Partial Differential Equations, Physiologically structured Models, Cell division cycle, Adaptive dynamics, Cancer, Drug control PACS: 87.17.Aa, 87.17.Ee

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Modelling targets for anticancer drug control optimisation in physiologically structured cell population models

Frédérique Billy * , Jean Clairambault * , Olivier Fercoq † , Tommaso Lorenzi * * , Alexander Lorz * and Benoît Perthame *

DRUG TARGETS IN PROLIFERATING CELL POPULATION DYNAMIC MODELS

McKendrick physiological age-structured cell population dynamics models for the cell division cycle [START_REF] Mckendrick | Proc. Edinburgh Math. Soc[END_REF][START_REF] Perthame | Transport Equations in Biology[END_REF][START_REF] Clairambault | A model of the cell cycle and its circadian control[END_REF] describe the evolution of densities n i (t, x) of cells of age x at time t in phase i of the cell cycle. Furthermore, one may include in such models targets for their physiological and pharmacological control [START_REF] Billy | [END_REF], either on death rates rate d i or on phase transition rates K i→i+1 (possibly set to zero by effects of cytotoxic drugs). A typical cell cycle model reads

         ∂ n i (t,x) ∂t + ∂ {v i (x)n i (t,x)} ∂ x + d i (t, x)n i (t, x) + K i→i+1 (t, x)n i (t, x) = 0, n i+1 (t, 0) = ∞ 0 K i→i+1 (t, x)n i (t, x)dx, n 1 (t, 0) = 2 ∞ 0 K I→1 (t, x)n I (t, x)dx. (1)
In a totally different setting, in an evolutionary dynamics representation, designed to describe the evolution of proliferating cell populations structured according to an evolutionary phenotype x, the dynamics of healthy (n H ) and tumour (n C ) cell populations undergoing a cytotoxic chemotherapy and subject to mutations may be written as

∂ ∂t n H (x,t)= growth with homeostasis 1 -θ H 1 + ρ(t) β r(x) - natural apoptosis d(x) - effect of drug c(t)µ H (x) n H (x,t) + θ H 1+ρ(t) β r(y)M σ H (y, x)n H (y,t)dy birth with mutation , (2)
∂ ∂t n C (x,t) = (1 -θ C ) r(x)-d(x)-c(t)µ C (x) n C (x,t) + θ C r(y)M σ C (y, x)n C (y,t)dy, (3)
where the total population is defined as

ρ(t) = ρ H (t) + ρ C (t), ρ H (t) = ∞ x=0 n H (x,t)dx, ρ C (t) = ∞ x=0 n C (x,t)dx. (4)
In this setting, drug effects have been put on the death term µ(x) only (this holds for cytototoxic drugs). But as will be presented below, it is possible to include the representation of cytostatic drugs as acting on the proliferation rate r(x).

DRUG EFFECTS IN HEALTHY AND CANCER CELL POPULATIONS Modelling drug effects on the cell division cycle, healthy and cancer cell populations

Solutions of linear models, such as McKendrick's, exhibit exponential growth. In Eq. ( 1), with time-periodic or stationary coefficients, solutions satisfy n i (t, x) ∼ C 0 N i (t, x)e λt [START_REF] Perthame | Transport Equations in Biology[END_REF][START_REF] Billy | [END_REF], if N i are eigenvectors of the eigenvalue problem

               ∂ N i (t,x) ∂t + ∂ {v i (x)N i (t,x)} ∂ x + λ N i (t, x) + d i (t, x)N i (t, x) + K i→i+1 (t, x)N i (t, x) = 0, N i+1 (t, 0) = ∞ 0 K i→i+1 (t, x)N i (t, x)dx, N 1 (t, 0) = 2 ∞ 0 K I→1 (t, x)N I (t, x)dx, N i > 0, N i (t + T, .) = N i (t, .), ∑ i T 0 ∞ 0 N i (t, x)dxdt = 1.
(

) 5
with T -periodic coefficients to account for physiological (circadian) and pharmacological (periodically delivered) control. The study of the growth exponent λ , first eigenvalue of this positive growth system, is therefore crucial. In [START_REF] Billy | [END_REF], we have shown how, taking transition rates K i→i+1 (t, x) as drug targets for a problem of optimal control by cytotoxic drugs on the eigenvalues λ H and λ C , population growth exponents in healthy and cancer cell populations, one may numerically solve the optimal control problem, using a Uzawa-like algorithm. This is illustrated on Fig. 1 and2. Moreover, possible targets for cytostatic drugs in this cell cycle setting are discussed in [START_REF] Billy | Designing proliferating cell population models with functional targets for control by anti-cancer drugs[END_REF]; they include velocities v i (x) (otherwise usually set to 1), the boundary term in the first phase (G 1 ) and the possible adjunction of a quiescent siding phase to which cells may escape at mitosis under the influence of a cytostatic drug, as in [START_REF] Gabriel | The contribution of age structure to cell population responses to targeted therapeutics[END_REF].

Modelling drug-induced drug resistance in cancer cell populations

In a setting slightly different from the general one presented above about evolutionary dynamics, we represent interactions of healthy and cancer cells separately with an environment variable. We assume that the cancer cell population is structured according to a 2-dimensional phenotype accounting for resistance to two categories of drugs: cytotoxic, acting on death rates, and cytostatic, acting on proliferation rates. Figures 3 and4 show that cancer cell populations, that in the absence of drugs develop a non-resistant phenotype more favourable to proliferation, may in the presence of the two drugs show evolution towards two different resistant phenotypes simultaneously [START_REF] Billy | Designing proliferating cell population models with functional targets for control by anti-cancer drugs[END_REF]. In [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], we show how this setting, without necessary inclusion of mutations, allows to set the bases for a representation of combined cytostatic and cytotoxic anti-cancer therapies that overcomes the development of drug resistances, and may lead to extinction of the cancer cell population while preserving the healthy cell population.

              
Within a general perspective of multidrug multiagent pharmacotherapeutic optimisation in oncology developed elsewhere [8,9], we have focused here on two different types of proliferating cell population dynamics models with physiological structure, both with possible inclusion of targets for cytotoxic and for cytotoxic drugs separately. In the first one, we have shown how the representation of cytotoxic drug effects simultaneously in healthy and cancer cell populations leads to an optimal control problem with unwanted toxic side effects constraints, which is numerically solved in [START_REF] Billy | [END_REF]; we note that with an accurate choice of targets for cytostatic drugs, it is possible to use the same numerical formalism to solve the same optimal control problem, using a combined cytotoxic-cytostatic therapy, as is frequently done in the clinic of cancers. In the second done, that belongs to the category of evolutionary dynamics models, already widespread in ecology but only recently developed in cancer biology [10], we present an adapted frame to theoretically tackle the problem of drug-induced drug resistance in cancer cell populations. We stress that in [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF] is shown a first numerical attempt at representing a combined cytotoxic-cytostatic therapy, that not only leads to avoid drug resistance, but even more, yields extinction of the cancer cell population while preserving the healthy one from irreversible damage.

We have focused here on cytotoxic and cytostatic drugs, i.e., on treatments targeting cell death and proliferation. It is noteworthy that a third fate awaits proliferating cells, which is differentiation, that is, maturation towards fully functional state. A characteristic of cancer cells is their insufficiently differentiated status, which may be due to a blockade of differentiation at some stage of their maturation process, as is the case in acute myeloblastic leukaemia (AML). It is known in the case of acute promyelocytic leukaemia (APL, or AML3, a variety of AML in which differentiation is blocked at the promyelocyte stage) that a combination of cytotoxic drugs and redifferentiating agents cures the disease in 80% of patients. A rationale for the use of such a combined treatment with targets for the two types of drugs is presented in [11], following a previous model of AML [12]. We hypothesize, following this case, that designing models adapted to specific cancer diseases and targeting the three main fates of cancer cells, i.e., proliferation, apoptosis and differentiation, and solving optimal control problems under constraints of toxicity and/or drug resistance, may in the future help clinicians to propose successful treatments in cancers that hitherto remain out of reach of standard chemotherapies in the clinic.

FIGURE 1 .

 1 FIGURE 1. Periodic (circadian) phase transition gating rates K i→i+1 (t, x) at G 1 /S (red dash-pointed lines), at G 2 /M without drug (green full lines) and at G 2 /M with multiplication by a cytotoxic drug effect (dashed blue line). Left: healthy cells; right: cancer cells. The optimal drug infusion scheme consists in rectifying the 'lazy' cancer gating into a stiffer one, close to the healthy gating. The cytotoxic drug is assumed to be S phase-specific, so that only the G 2 /M gating is affected by the drug. Redrawn from[START_REF] Billy | [END_REF].

FIGURE 2 .

 2 FIGURE 2. Effect of the drug on gating at G 2 /M transition. Left: without drug, the cancer cell population (blue line) outcompetes the healthy cell one (green line). Right: with partial drug blockade at G 2 /M, both populations are hit, but the healthy cell population wins over the cancer cell one. Redrawn from [4].

FIGURE 3 .

 3 FIGURE 3. Cancer cells competing with healthy cells in the absence of drugs: evolution of the cancer cell population from an average phenotype (x, y) = (0.5, 0.5) towards sensitivity to both drugs, the cytostatic (c κσ ) and the cytotoxic (c κτ ).

FIGURE 4 .

 4 FIGURE 4. Cancer cells competing with healthy cells in the presence of both cytotoxic and cytostatic drugs: evolution of the cancer cell population towards two resistant subclones, one resistant to the cytotoxic, the other one resistant to the cytostatic drug.