Michael Lange
email: michael.lange@imperial.ac.uk

Tony Field
email: t.field@imperial.ac.uk

Accelerating Agent-Based Ecosystem Models using the Cell Broadband Engine

This paper investigates how the parallel streaming capabilities of the Cell Broadband Engine can be used to speed up a class of agent-based plankton models generated from a domain-specific model compiler called the Virtual Ecology Workbench (VEW). We show that excellent speed-ups over a conventional x86 platform can be achieved for the agent update loop. We also show that scalability of the application as a whole is limited by the need to perform particle management, which splits and merges agents in order to keep the global agent count within specified bounds. Furthermore, we identify the size of the PPE L2 cache as the main hardware limitation for this process and give an indication of how to perform the required searches more efficiently.

Introduction

Multi-core stream processors have been used to provide significant performance gains for scientific algorithms in recent years [START_REF] Erez | Executing irregular scientific applications on stream architectures[END_REF]. The Cell Broadband Engine in particular has shown great promise for accelerating scientific applications due to its Stream Processing architecture [START_REF] Spray | Performance of a Lattice Quantum Chronodynamics kernel on the Cell processor[END_REF] [START_REF] Williams | The Potential of the Cell Processor for Scientific Computing[END_REF]. This paper investigates the use of the Cell Broadband Engine to accelerate a class of agent-based Virtual Ecology models for use in oceanographic research.

We focus on models generated by the Virtual Ecology Workbench, a domainspecific model compiler used by oceanographers to build and analyze models of plankton ecosystems in the upper ocean. VEW models have been used to study complex cause-effect relationships in ocean ecosystems, contributing to our understanding of the biodiversity of the ocean. The effect of marine plankton on the atmosphere and global climate, as well as the plankton ecosystem's response to external influences, like fishing or pollution, have also been the focus of VEWassisted studies [START_REF] Sinerchia | Testing theories on fisheries recruitment[END_REF] [START_REF] Hinsley | Planktonica: A system for doing biological oceanography by computer[END_REF].

Models generated by the VEW simulate large numbers of plankton particles following the Lagrangian Ensemble (LE) metamodel [START_REF] Woods | The Lagrangian Ensemble metamodel for simulating plankton ecosystems[END_REF], an individual-based approach that uses agents to simulate the biological and bio-chemical behaviour of individual plankters. In contrast to traditional population-based models this allows for a detailed analysis of life-histories of individual micro-organisms. LE modelling further facilitates the study of emergent properties of the ecosystem as a whole and avoids chaotic instabilities often observed in population-based models [START_REF] Woods | Stability and predictability of a virtual plankton ecosystem created with an individual-based model[END_REF].

The main workload of VEW-generated models involves updating large numbers of agents during each simulation timestep. Since the updates are independent of each other the process can, in principle, benefit greatly from parallel processing. Attempts to accelerate VEW models using FPGAs [START_REF] Lamoureux | Accelerating a Virtual Ecology Model with FPGAs[END_REF] have shown significant speedups, prompting further investigation of dedicated hardware solutions for parallel VEW simulations.

This paper investigates the use of Stream Processing principles to accelerate VEW models on the Cell processor, as well as the performance limitations imposed by agent management functions used in the VEW algorithm. We focus on efficient parallel agent update computation and make the following contributions:

-We describe a framework for parallel agent update on the Cell's SPE accelerator cores in Section 4. Our design includes buffered processing of agent data streams, as well as meta-data exchange to provide centralized scheduling. In addition to an implemented prototype simulation we also outline an extension to the VEW model compiler that generates SPU-specific SIMD vector code for performing agent update arithmetic. -An evaluation of the achieved performance gains based on a simple characteristic VEW simulation is given in Section 5. We show that the parallel speedups observed for the agent update loop, although promising, are limited by sequential overheads due to agent management functions. -We present an analysis of the observed limitations in Section 6. This shows that bottlenecks are caused by inefficient cache utilization. In this section we also present a hybrid processing approach that overcomes the identified limitations by exchanging meta-data between SPE and PPE processors to handle data searches and copy operations more efficiently. We demonstrate this method on a simple sequential search routine imposed through one of our design choices and largely eliminate the sequential overhead. -We discuss the implications of our findings for future research in Section 7.

In particular we focus on the current particle management function which is shown to limit parallel scalability on platforms such as the Cell.

Background

The Virtual Ecology Workbench (VEW) is a software tool designed to aid the creation and analysis of virtual plankton ecosystems based on the Lagrangian Ensemble (LE) metamodel [START_REF] Woods | The Lagrangian Ensemble metamodel for simulating plankton ecosystems[END_REF]. These models are individual-based and use agents (sometimes referred to as 'particles') to simulate large numbers of plankters and the various feedback processes needed to study oceanic ecosystems and their various emergent demographic properties. Since it is impossible to model every organism in the ocean individually, each agent in LE models represents a sub-population of identical plankton particles modelled after the behaviour of one individual plankter. Each agent follows its own trajectory and keeps an associated sub-population size in addition to internal state variables that are application specific. The sub-population size is used to infer demographic properties of the ecosystem, such as the global population count and the concentration of a particular species in a given domain.

The VEW provides oceanographers with an easy way of specifying classes of plankton species in terms of primitive phenotypic equations which are scientifically sound in that they are based on reproducible laboratory experiments. These equations define the basic behaviour for individual plankters in response to local environmental properties [START_REF] Hinsley | Creating Individual-based Models of the Plankton Ecosystem[END_REF], for example light, temperature, nutrient concentration etc. An agent's internal state is then updated using equations relevant to its current biological state.

After specifying plankton species and additional environment parameterizations, the VEW automatically creates code from a high-level domain-specific mathematical modelling language called Planktonica [START_REF] Hinsley | Planktonica: A system for doing biological oceanography by computer[END_REF] and provides further tools for the analysis of output data. The VEW framework substantially raises the level of abstraction in model development and facilitates the creation and analysis of complex virtual ecosystems without the need for conventional programming. It is an example of a domain-specific problem solving environment.

VEW models simulate virtual plankton ecosystems in a one-dimensional virtual water column that extends from the ocean surface to a specified depth (typically 500m). The column is divided into layers of 1m depth and currently ignores horizontal water fluxes; one imagines the column to be sealed on all sides with an open boundary at the bottom. The upper metre of the column is sub-divided into more fine-grained layers in order to accurately model light absorption near the sea surface.

The combination of solar heating and cooling to the atmosphere, combined with atmospheric wind stresses and other environmental processes determine the depth of the turbulent Mixing Layer in the upper ocean -the so-called Turbocline. The depth of the turbocline varies both diurnally and seasonally. Particles above the turbocline are subject to turbulent mixing, which is approximated in the model by random displacement1 . Below the turbocline there is laminar flow; here agents sink under the influence of gravity, although they may subsequently be 're-entrained' into the mixing layer, should the turbocline later deepen sufficiently to catch up with them. The dynamics of individual plankters is heavily influenced by turbulence.

Phytoplankton obtain energy from the sun through photosynthesis and absorb light in the process. Collectively, the plankton influence the profile of light and temperature in the column and this, in turn, affects the physics of turbulence -a process known as biofeedback.

VEW-created models simulate inter-agent processes, such as predation and ingestion, although these are not modelled using interactions between individual agents (an O(n 2 ) process). Instead, the population of a particular plankton type is aggregated into a field (a concentration for each layer); individual agents then interact with these fields (an O(n) process). This is the essence of the Lagrangian Ensemble approach.

Although agent updates can, in principle, proceed independently, the various fields (including nutrient fields) are shared. Any changes to these fields resulting from predation and the uptake or release of nutrients thus has to be negotiated among the agents. Furthermore, predation may reduce the concentration of a particular species which must, in turn, be reflected by a change in the agent's internal population counts. These processes are carefully managed in the generated code so that corrections from the previous time step and updates relevant to the current time step are merged into a single piece of 'update' code for each agent type.

The demographic noise, i.e. the statistical variability in the populations observed among identical runs of the model but with different random numbers, is limited by ensuring that a sufficient number of agents exists in the water column. For this purpose a Particle Management (PM) function is executed between timesteps. This splits agents with the largest sub-populations, creating new agents with independent trajectory. Similarly, in order to limit the computational cost of an over-populated simulation, the PM may also merge the least populated agents. The PM process thus provides a trade-off between computational cost and statistical accuracy for agent-based models.

Hardware

The Cell Broadband Engine is a heterogeneous stream processor which features a unique cache hierarchy. It consists of a central general-purpose CPU, the Power Processing Element (PPE), and eight vector processing nodes called Synergistic Processing Elements (SPE). These are much simpler in their architecture, but are optimized for high-throughput floating point arithmetic.

The PPE is a dual-threaded single-core general-purpose processor based on IBM's PowerPC architecture. It uses a typical CPU cache hierarchy including a 512KB L2 write-back cache with 128-bit cache lines. In our simulation the PPE will delegate the main workload of updating agents to the SPE arithmetic cores, whilst maintaining attributes of the water column in memory. In addition it is also responsible for the particle management task of the algorithm.

A circular on-chip data bus connects the PPE with the SPEs and transports data streams from main memory to the SPE accelerator nodes. The bus provides the high data bandwidth required for Stream Processing and has been shown to transport up to 196GB/sec in experiments with near perfect utilization [START_REF] Chen | Cell Broadband Engine Architecture and its first implementation[END_REF]. In practice, however, memory bandwidth is limited by the 25.6GB/sec XDR main memory interface.

The SPE accelerator nodes provide the computational power of the Cell through vectorized SIMD computation at 3.2GHz core frequency. Each SPE comprises a 256KB Local Store (LS) memory area for buffering stream data, as well as executable code. A separate Memory Flow Controller (MFC) coordinates asynchronous data exchange with main memory for efficient data localization.

The Synergistic Processing Unit (SPU) is the functional core of the SPE. It is optimized for SIMD vector arithmetic on vectors of four single-precision floats and uses instruction pre-fetching on two separate pipelines for Load/Store and arithmetic instructions. Although it is possible to use double-precision floats, these are not implemented natively and are therefore more costly. The original VEW works with double-precision, but we decided to only use single-precision in our prototype in order to investigate the performance limitations of the Cell processor available in the Playstation3.

The SPU uses branch-prediction in order to optimize loop constructs. The compiler predicts continued execution of the loop body, resulting in a 19 cycle penalty to exit the loop due to a pipeline flush. In order to avoid unnecessary pipeline flushes for conditional vector statements, outcomes for both branches are typically pre-computed. The results are then interleaved according to a vector of flags indicating the outcome of the condition for each element.

This particular style of execution demands a special-purpose instruction set, provided by IBM's Cell SDK through intrinsic functions. These include memory flow controls for asynchronous Direct Memory Access (DMA) via the MFC, message-passing primitives for communication between PPE and SPE (mailboxes), and vector instructions for arithmetic calculations and branch prediction.

Implementation

The main body of a VEW simulation comprises an agent update loop, which is the main source of parallelism in the model, and various housekeeping tasks that are predominantly sequential. Processing agent updates in parallel requires an efficient means of localizing agent data to the local SPE memory, as well as the exploitation of SIMD vector code for performing agent updates.

Memory Layout

Two important attributes of each agent are its associated species and its growth stage; these collectively define its type. There is a different update kernel for each type. To prepare agents for parallel update on the Cell they are stored as type-homogeneous agent arrays in memory.

Agent update is complicated by the fact that an agent may change its type at the end of a time step. For example, it may evolve to the next growth stage or it might die. Furthermore, agents may spawn new agents, for example as a result of reproduction. These various processes mean that agents may have to be moved to a different array before the start of the next timestep. This requires an additional search through each agent array during the sequential part of the algorithm. However, we will use this additional search routine to demonstrate how to combine the capabilities of PPE and SPEs to solve agent management tasks more elegantly (Section 6).

Due to the Cell's 128-bit bus lanes agents need to be aligned to multiples of 4 floating point numbers. In the prototype simulation an agent's state consists of 13 variables, requiring 16 floats to be allocated and stored in main memory. Conceptually this corresponds to an Array of Structures (AOS) layout (Figure 1 Vector processing on the SPU, on the other hand, requires a Structure of Arrays (SOA) layout, where each vector holds the same variable from four respective agents (Figure 1(b)). After buffering agents in the LS we therefore need to convert the agent state variables, as well as the environment variables, to vector float types in a small loop.

Since we are mostly processing four agents at a time, we can unroll this conversion loop, which allows the compiler to optimize the process through inlining. For conversion of less than four agents the loop accesses each element in all vectors in turn, allowing us to ignore the unused elements in incomplete vectors. This scheme is similar to a pad-with-zeros approach and is also used when preparing the data for export to main memory.

Agent Allocation

Since our framework aims to utilize the SPEs as pure Stream Processing nodes, the PPE handles all agent allocation and the scheduling of update tasks to the SPEs. This also provides the possibility of investigating different load-balancing schemes. A descriptor data structure is used to transport meta-data to the SPEs, which contains the memory address and size of each agent block. A twodimensional array of descriptor meta-objects is held in memory, mirroring the actual maintained agent blocks one-to-one. In addition we can use similar feedback objects to send meta-data back to main memory after processing a block of agents. This meta-data exchange allows us to control SPE processing in a non-static fashion and write an independent scheduling function that allocates a descriptor block to each SPE. We found an Round-Robin allocation scheme for independent agent arrays to provide sufficient load-balance between SPE processing nodes.

Agent Localization

The implemented data transfer framework follows general Stream Processing principles by masking all data transfer latencies through asynchronous DMA. By using three agent buffers we allow GET and PUT transfers to happen in parallel to SPU computation, allowing for continuous SPU execution.

This Triple-Buffering approach proved sufficient to handle all data transfers without incurring unexpected data stalls. The actual size of agent blocks showed no noticeable effect on DMA transfer times, as long as a minimum of 2 agents per block is maintained.

Since Stream Processing requires the agent data to be pre-fetched, we need to load the corresponding descriptor block before initializing the data GET. Thus we use a 2-step-lookahead loop to first fetch the meta-information, followed by the agent data, whilst maintaining the triple-buffered loop structure. This is illustrated in Figure 2. The different data and meta buffers are synchronized over a common counter (n).

Agent Updates

We extended the existing VEW code generator module to produce optimized SPU agent update code that can be compiled and executed with the designed agent transfer framework. This allowed us to generate inlined vector intrinsics resulting in fast agent update kernels with low-level optimizations.

Particular care needs to be taken when generating conditional statements, where we used pre-calculation of both branches, as described in Section 3. Through recursive compilation the model compiler is capable of handling nested conditionals of any depth. The non-standard architecture of the Cell makes it hard to compare accurately its performance to conventional x86 CPU architectures. For our performance evaluation we therefore used a 1.6GHz dual-core, as well as a 3.2GHz quad-core x86 CPU to run a single-threaded sequential version of our prototype simulation.

Performance Evaluation

The Cell implementation was developed and tested on a Playstation3 with 6 enabled SPEs, running Ubuntu 8.10 with IBM's Cell SDK and compilers2 . A single-precision floating point version of the prototype model was used for both architectures 3 .

For the purpose of this investigation we used a simplified version of the LERM [START_REF] Sinerchia | Testing theories on fisheries recruitment[END_REF] model, which simulates a Diatom (a type of marine phytoplankton) population in the Azores region for 2 years. This so-called Toymodel simulation incorporates all major components of the VEW algorithm and was handcoded for the Cell. For comparison we used an identical single-threaded sequential version written in C on the x86 platforms.

We scaled the size of the simulation by increasing the initial number of agents. Considering just the agent updates (i.e. turning off particle management), we observed average parallel speedups of more than 22 over the 1.6GHz runs of the simulation and more than 6 over the 3.2GHz version (Table 1). Furthermore, we verified that agent update runtimes scale approximately linearly with increasing numbers of SPE cores.

Figure 3(a) shows the total execution times, i.e. including particle management and various other sequential components. The graphs show the Cell is significantly faster than the 1.6GHz CPU but fails to provide significant overall speedup over the 3.2GHz processor. The overall non-linear scalability of the Cell version of the code is also quite evident. This is due to the particle manager, which is the most significant sequential component running on the PPE, as shown in Figure 3

More on Particle Management

Particle management consists of two main search functions, Split and Merge, that exhaustively search the agent arrays in order to identify the agents that need to be altered. The Merge function furthermore includes a O(n 2 ) lookup in order to find the relevant pairs for merging. The overall runtime shows a non-linear trend, although there is a linear number invocations of the Split and Merge functions.

Due to the AOS layout in main memory copying agents can be handled very efficiently on the PPE, since the L2 cache will already contain the whole agent once the first field is accessed. Similarly, the copied agent will be written to memory in one block due to the write-back property of the cache. On the other hand, searching through all agents whilst inspecting only one variable will be very inefficient since each agent is still loaded anyway. However, once the simulation contains more agents than the L2 cache can hold, each agent will be read multiple times during the Merge process(O(n 2 )), resulting in poor cache performance.

The size of the PPE's L2 cache therefore poses the main hardware limitation. The cache can hold 512KB which corresponds to 8K agents of 16 floats (64B), in contrast to the 2MB cache found on the x86 CPU, which may hold up to 32K agents.

In the x86 version of the code each timestep is dominated by the agent update loop, so the overhead of the particle manager is less evident. In the Cell version the cost of the update loop is substantially reduced and the cost of sequential computations on the PPE is increased when compared with an x86, as described above. The two combined have the effect of increasing the dominance of particle manager (see also Figure 4(b)), as per Amdahl's law.

Agent Type Changes

The same cache limitations can be observed when handling agent type changes on the PPE. As described in Section 4, an agent may arbitrarily change its growth stage during the update process. Thus we need to move all agents that did change type to the correct sub-array to ensure type-homogeneity before the next timestep. In contrast to particle management, however, we only need to traverse each sub-array once to identify all agents to be moved. We can therefore flag all state changes on the SPEs as they happen, and use the feedback metadata objects to send back the indices of the affected agents. Knowing the exact location of all agents that need to be moved, we can then utilize the PPE cache more efficiently to copy the agents, avoiding irregular memory access patterns. The effect of this method can be seen from Figure 4, which shows the respective contribution of parallel and sequential components. Figure 4(a) highlights the significant amount of time spent performing the search on the PPE, while Figure 4(b) shows improvements resulting from the parallel search strategy on the SPEs. This demonstrates that the SPE-assisted agent copy method substantially reduces the overhead imposed by this component, indicating how we can utilize the advantages of both types of processing core to perform agent management tasks very efficiently.

Unfortunately, the same ideas can not be used to optimize the particle management code. This is because the current management rules limit the total agent count globally. During parallel processing each SPE node only sees a local subset of all agents, and is therefore not able to tell how many agents it has to flag for Split/Merge copy operations.

Future Work

The original intention of this work was to explore ways to parallelize VEW models with hope of providing linear scalability with increasing numbers of agents. This was done in in order to prepare the VEW algorithm to be applied to three dimensional scenarios, which will inevitably feature agent counts several orders of magnitude larger than the existing models.

The greatest obstacle to parallelizing the VEW algorithm is the particle manager. The current rules are applied globally and thus do not scale well. We are currently exploring alternative management strategies that predominantly involve local computations and which avoid the need for sophisticated search strategies, e.g. to find the k largest or smallest agents by population count. The evaluation of new particle management algorithms involves an extensive combination of performance analysis and statistical analysis that is by no means trivial.

Conclusion

We have shown how to utilize Stream Processing capabilities of the Cell Broadband Engine to parallelize agent-based ecosystem models based upon the Lagrangian Ensemble metamodel. We have demonstrated parallel speedups for agent update computation of over 22 and 6 when compared to a 1.6GHz and 3.2GHz reference CPU respectively. We have also shown that the Cell can handle models with significantly large numbers of agents. The performance gains achieved are similar to speedups achieved for other scientific codes [8] [4] and indicate that the Cell processor has considerable potential for accelerating this class of applications.

We have found that the Cell DMA system can be used effectively to perform meta-data exchange in addition to agent data streaming. This adds flexibility to the overall design and allows us to perform SPE-based searches in conjunction with agent update computation. This approach overcomes one of the major drawbacks of SPE computation, which has been shown to work inefficiently for SPE-based searches without sufficient additional workload [START_REF] Bader | On the Design and Analysis of Irregular Algorithms on the Cell Processor: A Case Study of List Ranking[END_REF].

Through our investigation we also highlight that the particle management process, which is currently implemented sequentially, becomes a significant bottleneck. This highlights the importance of parallelizing the particle management algorithm.

Our investigation also demonstrates that the size of the PPE L2 cache can significantly limit performance. Since scientific applications often require several different memory access patterns the heterogeneity of the Cell's individual cores can be a great advantage. For example, randomly copying coherent data structures such as LE agents profits from the conventional cache hierarchy of the PPE, while exhaustive searches can be performed efficiently on the SPEs if combined with sufficient arithmetic workload. However, as shown in the previous section, the current implementation of the PPE does not compare well to conventional general-purpose CPUs, limiting its use as an additional processing core.

  (a)), which is maintained during the DMA transfer to the SPEs.

Fig. 1 .

 1 Fig. 1. Agent data layout in global and local memory.

Fig. 3 .

 3 Fig. 3. Execution times of parallel and sequential prototype.

  (b).

Fig. 4 .

 4 Fig. 4. Percentage of runtime spent in sequential and parallel parts of the main loop.

Table 1 .

 1 Runtimes of agent updates on Cell and x86 with respective parallel speedup.

	Agents Parallel runtime Sequential runtime Parallel speedup
		Cell	1.6GHz 3.2GHz 1.6GHz 3.2GHz
	4000	5.4s	96.7s	30.5s	17.91	5.65
	8000	9.5s	187.6s	58s	19.75	6.11
	16000	17.5s	368.3s	114.3s	21.05	6.53
	32000	34s	752.2s	230.8s	22.12	6.79
	64000	66.7s	1478.7s 454.1s	22.17	6.81

This approximation is good provided the time step is sufficiently large; we use a 30-minute time step

A PowerXCell i8, which provides improved double-precision floating point performance, was not available to us during this project.

[START_REF] Erez | Executing irregular scientific applications on stream architectures[END_REF] The use of single-precision only adds a small amount of random noise to an already noisy simulation.


Power Management Unit

Find supply (V i ) and threshold (V ti ) voltage assignments Such that application energy consumption is minimized Such that application energy consumption is minimized

( ) ( ) ( ) ∑ ∑ ∑ ∀ ∀ ∀ + = T i T i bit T i ti i i App j , i E j , i vol V , V E E min ∈ ∀ ∈ ∀ ∈ ∀ T i T i T i
Energy consumed when the task is executed at (V i V ti )

Communication energy consumption

Subject to the following deadline constraints per task t: This can be also implemented as a branch & bound algorithm.

x is executed at (V i ,V ti ) consumption
We can obtain exact results for small examples.

Voltage Assignment Algorithm

For all pairs of neighboring islands (i , j ) neighboring islands (i , j )

Merge VFIs i and j

Solve static VF assignment problem

Given an initial partitioning with N islands, find the static voltages

Compute the energy consumption

Merge the pair of islands that Update the VFI configuration

Merge the pair of islands that provides the minimum energy

( ) ( ) ( ) ∑ ∑ ∑ ∈ ∀ ∈ ∀ ∈ ∀ + = T i T i bit T i ti i i App j , i E j , i vol V , V E E min x RGM2-ISCA
29% 0 5 J 2-VFI: 7.5mJ 3-VFI: 7.6mJ 0 1 2 0 1 2 0 Su ( , ) tile (2,1) tile (1,1) tile (1,2) Single VFI 29% 1.5 V) 0 0 Single VFI 1.5 V) 0.5 1 ply Voltage (V 0.5 1 ply Voltage (V 1 2 1 2 0 Supp tile (1,1) tile (1,2) tile (2,1) tile (2,2) 1 2 1 2 0 Supp tile (1,1) tile (1,2) tile (2,1) tile (2,2) RGM2-ISCA'10 110 0 1 0 1 Two VFIs 0 1 0 Three VFIs

Experiments with Realistic Benchmarks Several E3S benchmarks (consumer, network, auto-industry, telecom)

A li i h d l d N C i f 3 3 Applications scheduled to NoCs ranging from 3×3 to 5×5

PE PE PE ] [ 2 1 N q , ... , q , q Q = State (queues utilization) PE PE PE ] [ 2 1 N f f f F = M Input (clock speeds) RGM2-ISCA'10 118 PE PE PE ] [ 2 1 N f , ... , f , f F M
Step-by-step Model Construction (one queue) Amount of data (packets) read from the queue

λ μ V 1 , f 1 λ λ μ V 1 , f 1 Average utilization in the k th control interval ) 1 ( ) 1 ( ) 1 ( ) ( 1 1 - - - + - = k Tμ k Tλ k q k q q 2 q 1 λ 1 μ 2 λ q 2 q 1 λ 1 λ 1 μ 2 λ
Amount of data (packets) written to the queue

μ 1 λ 2 V 2 , f 2 μ 1 μ 1 λ 2 V 2 , f 2
If data read/write rates are proportional to the frequency of the VFI

) 1 ( ) 1 ( , ) 1 ( ) 1 ( 2 1 1 1 1 1 - = - - = - k f μ k μ k f λ k λ ) ( ) ( , ) ( ) ( 2 1 1 1 1 1 f μ μ f
The state-space equation can be written as

[ ] ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - - + - = ) 1 ( ) 1 ( ) 1 ( ) ( 1 1 1 k f k f μ λ T k q k q RGM2-ISCA'10 119 ⎥ ⎦ ⎢ ⎣ -) 1 ( 2 k f B
Step-by-step Model Construction (three queues)

⎥ ⎤ ⎢ ⎡ -1 1 0 μ λ V 1 , f 1 V 2 , f 2 λ 1 μ First row → q 1 ⎥ ⎥ ⎥ ⎥ ⎢ ⎢ ⎢ ⎢ - = 2 2 0 λ μ B q 1 q 2 q 3 1 μ 1 μ 2 λ 3 Second row → q 2 ⎥ ⎦ ⎢ ⎣ -3 3 0 μ λ V 3 , f 3 λ 2 μ 3 f f f Third row → q 3
The topology of the VFIs determines the matrix B

f 1 f 2 f 3

An algorithm automatically constructs B

The structure of the model is the same regardless of B

1 1 1 ) 1 ( ) 1 ( ) ( × × × × + = M M N N N k- F TB k- Q k Q RGM2-ISCA'10 120 1 1 1 × × × × M M N N N

System Controllability

In the multiple voltage-frequency island system with M islands, utilization of at most M queues can be controlled. The system is controllable iff rank(B) = N (i.e., number of controlled queues)

Motion

Feedback Control Architecture

R(k) K

∑ + B ∑ + 1 I Q(k) Gain matrix F (k) R(k) K 0 ∑ B ∑ z -1 I N×N Q(k) Desired Queue Utilizations + I N×N Utilizations Queues under control ) 1 ( ) 1 ( ) ( k TBF k Q k Q + = Open loop system: K State feedback matrix ) 1 ( ) 1 ( ) ( k- TBF k- Q k Q + = State feedback matrix Cl d l t ⎥ ⎤ ⎢ ⎡ -1 1 0 μ λ 0 ) 1 ( ) ( ) ( TBRK k- Q TBK I k Q + - = Closed loop system: ⎥ ⎥ ⎥ ⎥ ⎢ ⎢ ⎢ ⎢ - = 2 2 0 0 λ λ μ B RGM2-ISCA'10 122 0 ) ( ) ( ) ( Q Q ⎥ ⎦ ⎢ ⎣ -3 3 0 μ λ

  i.e., optimum number of VFIs, n ≤ N) Assignment of the supply and threshold voltages to each island Assignment of the supply and threshold voltages to each island Such that the total energy consumption is minimized

Finds

Partitioning and Voltage Assignment Algorithm For all pairs of neighboring islands (i j ) neighboring islands (i , j ) Merge VFIs i and j Solve static VF assignment problem Merge VFIs i and j Given an initial partitioning with N islands, find the assignment problem Compute the ti with N islands, find the static voltages energy consumption Update the VFI Merge the pair of islands that Solve the static voltage/frequency assignment problem Update the VFI configuration g p provides the minimum energy assignment problem RGM2-ISCA'10 107

	VFI
	f x	t t	+	t	t Comm	≤	deadline	t	-	start	-	time	t
	Execution time			Communication delay
	RGM2-ISCA'10												106

'10 108 t t t Comm t t time start deadline t f x

	subject	+	≤	-	-

to Voltage Assignment Algorithm For all pairs of neighboring islands (i j ) neighboring islands (i , j ) Merge VFIs i and j Solve static VF assignment problem Given an initial partitioning with N islands, find the static voltages Compute the energy consumption Merge the pair of islands that Update the VFI configuration Merge the pair of islands that provides the minimum energy Constrained nonlinear optimization or nonlinear programming

VFI solution Outline VFI partitioning Multi-VFI NoC designs P titi i d lt i t Partitioning and voltage assignment Examples On-line control State-based model construction Feedback control architecture Feedback control architecture Stability issues Summary Summary RGM2-ISCA'10 Why On-line Control? Cannot rely on nominal values because they vary Temperature because they vary Sources of concern are workload, process, voltage, temperature

	20 16 20 8 12 4 8 0 Improvement over tion (mJ) Consump al Energy Tota RGM2-ISCA'10 p variations consumer 36% g 1-Delay Minimum energy gy network auto-industry telecom 1 VFI 1-VFI 2-VFI 3-VFI 49% 80% 78% Target and actual queue ili i p Cope with	111

the parameter variations which cannot be Controller utilizations variations which cannot be predicted or accurately modeled at design time System under control Voltage Frequency Heuristic techniques and manual tuning won't work! Static Power Dynamic Power y RGM2-ISCA'10 113 Distributed Power Management in Magali A 477mW NoC-based digital baseband for MIMO 4G SDR chip organized RGM2-ISCA'10 114 around a 15-router asynchronous NoC that connects 22 processing units. [Clermidy et al, ISSCC'10] Local Control in Multi Clock Domain Processors The clock domain partitions in an MCD processor Interface model between the domains

	[Semeraro, et al, HPCA'02]	[Wu, et al, ASPLOS'04]

PID controller for voltage/frequency control proposed previously using only local queue information Ignores interactions among multiple queues Works fine if frequency change in one clock domain has negligible impact on other domains For an MCD processor with arbitrary partitions and strong interactions among multiple RGM2-ISCA'10 115 For an MCD processor with arbitrary partitions and strong interactions among multiple queues, a centralized online DVFS scheme may be needed Design Methodology for Multi-VFI NoCs Traditionally, PID controllers are used due to simplicity. However, state-space modeling brings new opportunities state space modeling brings new opportunities

	Precise controllability and stability analysis
	Pole placement linear quadratic regulator robust controller Pole placement, linear quadratic regulator, robust controller
	ns utilization Desired u D	Os face FIFO for interf	+ _	∑ ∑	Voltage-frequency q y controller	V 1 , f 1 V 1 , f 1 V 2 , f 2 … V N , f N …	NoC under control
						State	Actual utilization of interface FIFOs

RGM2-ISCA'

10 116 feedback of interface FIFOs Design Methodology for Multi-VFI NoCs Parameter Variations Workload VFI Configuration State-space Model Construction

							Feedback Control
	Configuration	Construction	Design Design
	lizations esired uti De	ce FIFOs or interfac fo	+ _	∑ ∑	Voltage-frequency f controller	V 1 , f 1 V 2 , f 2 V N , f N …	NoC under control
						State	Actual utilization of interface FIFOs
	RGM2-ISCA'

10 117 State feedback of interface FIFOs Formal Feedback Control Multi-VFI Network-on-Chip Write to FIFO (λ packets/sec) Operating frequency f 2 Interface queue utilizations are the states of the system State feedback for voltage (λ packets/sec) Read from FIFO Operating f f frequency f 2 State feedback for voltage- frequency control Control interval is T μsec

	frequency f 1	Read from FIFO (μ packets/sec)
	μ
	State (queues utilization)

Quantization Inv. DCT Buffer Frequency Tracking Capabilities 50 Frames/sec for 352×288 CIF frames f 1 is set to meet the target, f 2 and f 3 follow f 1

	y requency Clock F	GHz) (G	0.1 0.2 0 0	Nominal clock frequency (f 2 ) 20 40 60	Actual clock frequency 80 100
	equency Clock Fre C	z) (GH	0.02 0.04 0 -0.02 0 0	C t l I t Control Interval l 40 60 Nominal clock frequency (f 3 ) 20 20 40 60	46% power savings 100 Actual clock frequency 80 80 100
					Control Interval
	RGM2-ISCA'10		126

Results on FPGA prototyping Work on FPGA prototype using Virtex-II Pro FPGA from Xilinx Inter domain communication Inter-domain communication Delay Locked Loops (DLLs) used to generate individual clock signals Block RAM based mi ed clock FIFOs Block-RAM based mixed-clock FIFOs Voltage conversion not supported yet by Xilinx boards MPEG-2 encoder design divided into three VFIs

	S h Synchronous design utilizes 16966 LUTs d i tili 16966 LUT Design with three VFIs utilizes 19161 LUTs	13% overhead
	P Power consumption obtained using XPower t i b ti d i X P
	Without voltage scaling, power drops from 277W to 259W
	Consistent with simulations Consistent with simulations
	RGM2-ISCA'10	127

Clock Control Architecture 4 basic clocks Clock Control Clock Source Clock Divider Algorithm ROM PicoBlaze Microprocessor 17.5 15 Clock DLLs Four basic

	Clock Control Architecture
						Clock Control	4 basic clocks
			Clock Source		Algorithm ROM
						PicoBlaze	Clock
						Microprocessor Microprocessor	Divider Divider
			Clock DLLs
	Four basic Four basic clocks	12.5 MHz	15 MHz	17.5 MHz	20 MHz	Search Search Frequency
			Clock Control		Frequency Table ROM ROM
	Clock 1 Clock 1		Clock 2 Clock 2			Clock 3 Clock 3
	PE1		PE1			PE1
	MicroBlaze		MicroBlaze		MicroBlaze
	Microprocessor Microprocessor	FIFO 1	Microprocessor Microprocessor	FIFO 2	Microprocess Microprocess or
	RGM2-ISCA'10						129
						Microprocessor	Divider
	Four basic clocks	12.5 MHz	MHz	MHz	20 MHz	Search Search Frequency
			Clock Control		Frequency Table ROM ROM
	Clock 1 Clock 1		Clock 2 Clock 2			Clock 3 Clock 3
	PE1		PE1			PE1
	MicroBlaze		MicroBlaze		MicroBlaze
	Microprocessor Microprocessor	FIFO 1	Microprocessor Microprocessor	FIFO 2	Microprocess Microprocess or
	RGM2-ISCA'10						128

Clock Control Architecture 4 basic clocks Clock Control Clock Source Clock Divider Algorithm ROM PicoBlaze Microprocessor 17.5 15 Clock DLLs Four basic

						Microprocessor	Divider
	Four basic clocks	12.5 MHz	MHz	MHz	20 MHz	Search Search Frequency
			Clock Control		Frequency Table ROM ROM
	Clock 1		Clock 2	The voltage/	Clock 3

frequency control circuitry can PE1 MicroBlaze Microprocessor Clock 1 PE1 MicroBlaze Microprocessor Clock 2 PE1 MicroBlaze Microprocess Clock 3 The voltage/frequency control circuitry can be implemented in HW The control algorithm can be implemented

	Microprocessor	FIFO 1	i fi Microprocessor in firmware / OS / OS FIFO 2	Microprocess or
	RGM2-ISCA'10			130

Feedback Control Architecture

Feedback Control Architecture

R(k) K

Experiments with MPEG-2 Encoder

The encoder is divided into three VFI islands and mixed clock FIFOs are used at the interfaces FIFOs are used at the interfaces The frequency of Variable Length Encoder is set to achieve the desired encoding rate g


Motohiro Takayama

Ryuji Sakai

Parallelization Strategy for CELL TV

Consumer electronics devices are moving forward to utilize multi-core processors. We have developed a series of unique applications for TV based on Cell Broadband Engine TM (Cell/B.E.) 1 . This paper introduces such applications realized by the capability of the multi-core processor, and shares the strategy we took to exploit its potential.

Introduction

CELL TV TM [CT] is a real world application utilizing a multi-core processor in the consumer electronics industry. It has unique media processing applications such as simultaneous playback of 8 channels, and 8 channel recording named Time Shift Machine. These applications are realized by the processing power and high bandwidth of the heterogeneous multi-core processor Cell/B.E.

In order to provide users with these attractive features, we have created a framework for developing media processing applications, and a parallel programming model for exploiting the parallelism in modules of the applications. We combined the framework with the programming model in a nested fashion. At the coarse level, the modules of an application are scheduled by the framework. The framework manages the concurrency and keeps the resource criteria. To parallelize the module that requires more than single core, the programming model is used. It eases the parallelization and gains the scalability. We also vectorized the serial portion of the module by SIMD instruction to exploit data parallelism.

Applications

The computation power and high bandwidth of Cell/B.E. makes the multiple media processing work simultaneously. CELL TV plays 8 2 channels simultaneously and composite into one display, and records every content of 8 channels within the last 26 hours. With simultaneous playback of 8 channels, users can see all broadcasted channels at once. It replaces the channel zapping. Users just select the channel they would like to watch, instead of pressing the next/previous channel button continuously. Time Shift Machine is a feature that records whole contents of every channel within the last 26 hours. That throws the idea "forget to watch" away. Users can watch any contents of any channels In addition, users can watch compute intensive YouTube HD video while recording all channels. This example shows the application is well parallelized to supply both required processing power and high bandwidth.

Those applications are comprised of several media processing modules running on SPEs. Example of those modules are audio/video codec (MPEG-2, H.264/AVC, AAC, etc), stream demultiplexer, and image enhancement filters (Super Resolution [SR], noise reduction, etc). Those modules are connected to compose a single functionality of a complex application.

Fig. 1 shows the connection of modules to play a audio/video stream as an example. Though a module may require more processing power than a single core offers, application developers do not have to care about it. This modularity simplifies complex application development. Fig. 2 depicts how the simultaneous playback of 2 streams works. It is obvious that the multiple stream playback is a simple extension to the single stream playback case.

Parallelization Strategies

To make the applications run in real-time, we took a three-level parallelization strategy: inter/inner module parallelization, and SIMD. Fig. 3 describes the overview of the strategy.

The inter/inner module parallelization strategy is analogous to the separation of the OS and the thread library. In ordinary PC, the OS schedules the process to provide the users a view that applications are running concurrently. On the other hand, the thread library has the responsibility to make the application run in parallel, and provides the platform portability.

Inter Module Parallelization

AV Application framework [SM] has a responsibility for concurrency among the modules. The framework provides the resource reservation facility and real-time task scheduling. The connection of modules is defined in an XML file with resource attributes. The attributes specify such resources as how long the module uses the processing core, how many cores the module uses, how many memories are necessary for I/O, etc. The framework harmonizes the modules running concurrently on SPEs with satisfying real-time criteria. It is just like a real-time OS for SPEs.

Inner Module Parallelization and SIMD

Although the framework manages the concurrency and schedules resources appropriately, there still remains a challenge in parallelizing compute intensive modules to achieve the required performance.

We have created a programming model named Molatomium [MT] to ease the parallelization. Programming in Molatomium is composed of two part: Mol and Atom. Mol is a coordinate language that describes the parallelism among Atoms. The code written in Mol is compiled into portable byte code and executed by virtual machines residing in SPEs. Mol is suitable for writing a parallel algorithm in task-parallel manner. On the other hand, Atom is a platform native code written in C/C++ like a kernel in CUDA [CU] and OpenCL [CL]. It is a parallel execution unit of Molatomium, and is optimized with platform dependent facility like SIMD instructions.

Although media processing applications are well suitable for using SIMD type data parallelization, some algorithms like Super Resolution do not really fit into this category. Molatomium can be applied to such case because it can not only describe data parallelism in Atom but also express task parallelism in Mol.

The parallelization strategy for CELL TV proved the hierarchical approach works well on the Cell/B.E. architecture, which has a programmable memory hierarchy between the local storage and the main memory. It suggests that the hierarchy parallelization strategy would be effective in the coming NUCA (Non Uniform Cache Architecture) platform as well as the processor core count increases.

Though the modular design of the application framework provides the independency and increases the reusability, the task assignment attributes used by the scheduler are static. For instance, while the calculation time the module given is suitable for one content, other contents may require less computation time. We are working on improving Molatomium capability so that it can describe a whole application. Though it lacks the ability to assert real-time criteria so far, It could improve core utilization by its dynamic load balancing. In contrast to AV Application framework which is designed especially for the requirements in developing Cell/B.E. platform, Molatomium provides platform portability. In order to expand the applications created for CELL TV to other platforms, the whole application would be described by Molatomium.

Though the framework and the model work fine, we still tend to get lost in debugging bandwidth issue. That issue will increase in the near future as the many-core processors come to main stream. To solve this, the hardware and the development tool should provide the way to optimize the bandwidth consumption with ease.

Conclusion

CELL TV is our latest challenge to apply the multi-core processor to the consumer electronics device. In order to meet the required performance, we parallelized the applications using three-level parallelization strategy: AV application framework for concurrency, Molatomium for parallelizing the modules, and SIMD for exploiting the data parallelism.

  Fig. 1. Stream Decoding

  Fig. 3. Hierarchical Parallelization Strategies

Cell Broadband Engine and Cell/B.E. are trademarks of Sony Computer Entertainment, Inc., in the United States, other countries, or both and is used under license therefrom.

 2 The TV supports 8 channels because there are 8 channels for digital terrestrial television in Japan.


Nicolas Chenouard

Alexandre Dufour

Dr Jean-Christophe Olivo-Marin

Tracking algorithms chase down pathogens

Keywords: particle tracking, cell tracking, pathogens, shape, motility

Understanding sub-cellular dynamic processes governing pathogenic mechanisms is a necessary step towards the development of new drugs and strategies against infectious diseases. Sub-cellular pathogenic mechanisms, such as viral invasion processes involve highly dynamic nano-metric scale objects and rapid molecular interactions that require the study of individual particle paths. Single particle tracking methods allow to visualize and characterize the dynamics of biological objects and provide a straightforward and accurate means to understand sub-cellular processes. This review describes a number of particle tracking methods in time-lapse microscopy sequences and provides examples of using such techniques to investigate mechanisms of host-pathogen interactions.

Pathogenic mechanisms are inherently governed by highly dynamic processes both at the molecular and the cellular level. Understanding and interfering with these processes are necessary steps toward the development of new drugs and strategies against infectious diseases. The path towards this goal has been significantly enlightened over the past decades thanks to the numerous advances in fluorescent probes, labeling techniques and microscopy systems, which allow one to look directly at the biological system at the appropriate scale both in space and time. Recent fluorescent imaging protocols have thus allowed scientists to shed light on a large amount of biological phenomena both at the cellular and sub-cellular level, many of which had never been observed in other imaging modalities [START_REF] Vonesch | The colored revolution of bioimaging[END_REF].

The study of a dynamic process is intricately related to its spatio-temporal scale, both on the imaging and on the quantification side. More specifically, the distinction between cellular and molecular scales has lead to different imaging techniques and different approaches to analyze the resulting images. This review focuses on nano scales while a companion paper in a future issue will cover cellular scales.

Sub-cellular pathogenic mechanisms, such as viral invasion and spreading of prions involve highly dynamic nano-metric scale objects and fast interactions. Indeed, pathogens usually exhibit rapid motion within cells, especially when supported by an active transport mediated through F-actin, myosin or motors of the microtubule network [START_REF] Lehmann | Actin-and myosindriven movement of viruses along filopodia precedes their entry into cells[END_REF]. For example, viruses can spread from cell to cell on small protrusions either driven by actin polymerization or along thin connecting fibers [START_REF] Brandenburg | Virus trafficking -learning from single-virus tracking[END_REF], and tunneling nanotubes may also represent a means of cell-to-cell communication used by prions for intercellular spread [START_REF] Gousset | Prions hijack tunnelling nanotubes for intercellular spread[END_REF]. The characterization of these fast processes is of major interest to understand

Early works

An intuitive idea to solve the tracking problem is to decouple a detection stage from a linking procedure between the measured positions. When imaging cellular processes, the fluorescently labeled objects often appear as small bright spots superimposed to an uneven background. A spot is defined as a relatively small and compact object with an intensity that is both diffuse and higher than the immediate neighborhood. Detecting spots in a noisy background is a fairly standard image processing task and many efficient methods have been designed to obtain a set of point coordinates that correspond to the objects of interest [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF][START_REF] Thomann | Automatic fluorescent tag detection in 3D with super-resolution: application to the analysis of chromosome movement[END_REF][START_REF] Aguet | A maximum-likelihood formalism for subresolution axial localization of fluorescent nanoparticles[END_REF][START_REF] Ober | Localization accuracy in single-molecule microscopy[END_REF]. However in the context of living cellular systems, spot tracking is made difficult by the poor quality of images and the fact that the dynamics of spots can change over time. Moreover when the density of objects is high, the targets can aggregate temporarily which makes the appearance of spots change and their localization ambiguous.

Conventional multiple particle tracking techniques rely on a detection step of the spots followed by a nearest neighbor association (NNA) [START_REF] Apgar | Multiple-particle tracking measurements of heterogeneities in solutions of actin filaments and actin bundles[END_REF] or a constrained NNA [START_REF] Crocker | Methods of digital video microscopy for colloidal studies[END_REF] to build the trajectories from the set of detections. The principle is fairly intuitive: in an iterative way a trajectory is extended with the nearest detection in the next frame. For a comparative study of these techniques for tracking multiple particles in 2D images refer to [START_REF] Cheezum | Quantitative comparison of algorithms for tracking single fluorescent particles[END_REF]. These kind of methods produce satisfactory results when image sequences show a limited number of objects that are easily identifiable and move slowly. However in practice these conditions are rarely satisfied, especially when imaging fast processes in 3D. So as to keep the acquisition time as short as possible, the exposure time is generally set very short, resulting in low intensity spots and a very high acquisition noise. In this case the detection process may fail and the detection of noisy artifacts can occur due to the low signal-to-noise ratio.

In order to overcome the problem of corrupted detections, methods that do not use a frame by frame detection correspondence have been proposed to perform simultaneously detection and tracking [START_REF] Sage | Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics[END_REF]. A 2D sequence in time is considered as a 3D volume and the tracking problem is formulated as a global optimization process solved by dynamic programming. In a given track, the method looks for the best path in the spatiotemporal volume that links high intensity pixels and has minimum displacements. This method is able to extract a particle trajectory even when images are extremely noisy. However A similar approach has been proposed in [START_REF] Bonneau | Single quantum dot tracking based on perceptual grouping using minimal paths in a spatiotemporal volume[END_REF]. Here the detection of particles and the creation of trajectories are decoupled tasks but the association is modeled as a perceptual grouping problem in the spatiotemporal volume. The tracking model is defined such that each individual trajectory is a minimal path in a Riemannian metric. The proposed algorithm extracts trajectories sequentially thanks to a fast marching procedure. The method is robust to transient object disappearance and was successfully used to track blinking quantum dots on membranes in 2D images.

Both methods use a regularity constrain of the spatiotemporal volume to extract trajectories even with high levels of noise or temporary disappearance of the particles, but a number of additional assumptions are made that are not always biologically relevant. For example, a strong requirement is that particles should have a slow motion compared to the acquisition rate, so that they are nearly overlapping in two successive frames and produce smooth spatiotemporal curves in the 2D+t volume. This is not adapted to most 3D applications in which movements are fast compared to the acquisition rate, and limits the use of the methods to 2D applications. Another similar approach has been proposed in [START_REF] Racine | Visualization and quantification of vesicle trafficking on a three-dimensional cytoskeleton network in living cells[END_REF] where the 3D+t space is projected in a 2D space called kymograph, where trajectories are extracted as smooth curves. This method has proved efficient on data sets where movements are linear and slow.

As a general fact, one can say that methods based on the regularity of trajectories are difficult to adapt to most biological applications. They aim at extracting long and smooth trajectories while in reality trajectories are often short and very irregular due to fast changes in motion. Moreover, prior knowledge such as the level and the distribution of noise in the images, the lifetime of trajectories, statistics of appearance or disappearance of the particles or probability of motion changes are difficult to incorporate into such models. Also, as their formalism lacks generality and flexibility, these methods are difficult to transpose out of the dedicated context they were developed for. In the sequel, we focus on the Bayesian filtering approach, which is a more general framework based on statistical modeling, and was first developed in the radar surveillance and computer vision communities.

Bayesian multiple target tracking

The Bayesian approach to tracking is inspired from the estimation theory [START_REF] Bar-Shalom | Tracking and data association[END_REF] and was extensively used for years in fields such as radar surveillance [START_REF] Bar-Shalom | Multitarget/Multisensor Tracking[END_REF] and computer vision. This approach consists in building the most probable trajectories given the set of detections in the entire sequence. The objective score function to maximize is defined as the likelihood of the association between tracks and detections. Corresponding probabilities are generally decomposed following the Bayes rule, which allows one to estimate the likelihood of association in an iterative way. In practice, most methods perform the computation sequentially in the order of detection occurrence, which is suitable to real time tracking applications.

In most cases trajectories can be considered independent, so the association likelihood is computed as a combination of the individual track probabilities, which are computed separately by adopting a so-called state space formalism [START_REF] Bar-Shalom | Tracking and data association[END_REF][START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF]. The state of a target is a vector containing its coordinates at a given time. In practice the state is unknown and we have only access to a measurement whose relationship with the state is defined in the state space by the measurement equation. The accuracy of the measurement procedure is limited by the resolution of the acquisition device and the accuracy of the localization algorithm. So in order to take into account a possible lack of precision a random term is included in the measurement equation. The evolution equation models the dynamics of a target in the state space. The equation defines a relationship between the current and previous states. As biological targets do not show a deterministic behavior but instead exhibit a random motion such as diffusive Brownian motion, a random term, coined process noise, is introduced in the evolution equation. In order to get relevant probabilistic scores from the dynamics model, the evolution equation has to model accurately the true behavior of the target. Diffusive motion has attracted a lot of attention in the biophysical community and many subtypes have been described, each with different models of motion [START_REF] Qian | Single particle tracking analysis of diffusion and flow in two-dimensional systems[END_REF][START_REF] Saxton | Single-particle tracking: Applications to membrane dynamics[END_REF]. Brownian motion, however, which is a special case of diffusive motion, is the most used in the biological tracking community. Statistical physics has shown that the position of a Brownian particle after a macroscopic time lag follows a multidimensional Gaussian distribution centered on the previous position [START_REF] Nelson | Dynamical theories of Brownian motion[END_REF], therefore the state evolution probability is written as a Gaussian function. The covariance matrix quantifies the motility of the particle: it depends on the time lag and the diffusion coefficient, which is influenced by physical data such as the shape of the target, its weight, the temperature and the viscosity of the medium.

Biological objects often move along a preferred direction when they evolve in an active flow or when they interact with the cytoskeleton, like, for instance, cargos on microtubules that move in a linear fashion with a high speed. In this case, directed or curvilinear models of motion are accurate. For most biological targets though, the motion is not known a priori, and, moreover, it can change during the observation time lapse. Hence it was proposed in [START_REF] Genovesio | Interacting Multiple Model based method to track moving fluorescent biological spots[END_REF] to describe the target motion as an adaptive mixture of three of the most general kinds of movement: Brownian, directed and curvilinear motions. Based on the Interacting Multiple Model (IMM) filter [START_REF] Blom | The interacting multiple model algorithm for systems with Markovian switching coefficients[END_REF], the weights of the model mixture self adapt to the observed motion of the target. By doing so it is possible to track targets for which there is little or no knowledge about the dynamics. The parameters of the motion models are generally fixed based on prior knowledge about the target behavior [START_REF] Saxton | Single-particle tracking: Applications to membrane dynamics[END_REF]. However it is also possible to define appropriate motion models that include an online estimation technique of these parameters. A good overview of such techniques is given in [START_REF] Li | Survey of maneuvering target tracking: dynamic models[END_REF].

Thanks to the evolution and measurement equations, the probability of a given measurement knowing all past detections generated by a target can be computed.

Generally, linear equations with Gaussian noises are appropriate. In this case, Kalman filtering [START_REF] Bar-Shalom | Tracking and data association[END_REF] is the optimal estimation technique. Techniques such as the extended Kalman filter [START_REF] Bar-Shalom | Multitarget/Multisensor Tracking[END_REF] or particle filtering [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] are used in the non linear case to approximate probability distributions.

In the field of cellular imaging, the Bayesian tracking model was first used in [START_REF] Ngoc | Adaptive detection for moving biological objects in video microscopy sequences[END_REF] to track particles in 2D and extended to 3D in [START_REF] Genovesio | Tracking of multiple fluorescent biological objects in three dimensional video microscopy[END_REF]. The principle of this approach is to use an iterative decomposition of the measurements likelihood to extend a set of active tracks frame by frame. The associations are sequentially chosen following a maximum likelihood principle: the most likely association between a track and a detection is selected, then the associated measurement and track are removed from the association problem. The procedure is repeated until there are no more tracks or detections to associate. The measurements that remain unassigned are used to instantiate new tracks, while a track with no detection can be extended during some frames thanks to a prediction given by the kinetic filter. The latter possibility models the temporary disappearance of a particle which is especially important when processing noisy images or when tracking blinking objects such as quantum dots. This algorithm has been used successfully in 3D biological applications such as HIV-1 virus tracking [START_REF] Arhel | Quantitative four-dimensional tracking of cytoplasmic and nuclear HIV-1 complexes[END_REF] and Golgi units tracking in the Drosophila oocyte [27]. It is worth pointing out that 3D imaging facilitates the tracking as it reduces the occlusion problems that occur in 2D images. Particle filtering has been proposed in [START_REF] Smal | Bayesian tracking for fluorescence microscopic imaging[END_REF][START_REF]Particle filtering for multiple object tracking in molecular cell biology[END_REF] as an alternative resolution scheme for the Bayesian tracking problem in 2D images and used in [START_REF] Godinez | Tracking of virus particles in time-lapse fluorescence microscopy image sequences[END_REF] to track HIV-1 particles.

Particle filtering is a sampling technique used to model the uncertainty of the association between tracks and targets. The algorithm maintains several hypotheses of association called particles in order to implicitly delay the association decision to a latter time. The selection of association hypotheses is achieved by a probabilistic resampling of the particles at each frame. An up-to-date review of particle filtering techniques can be found in [START_REF] Johansen | A tutorial on particle filtering and smoothing: fifteen years later[END_REF]. There are however a number of critical issues like the selection of the optimal number of particles, the sample impoverishment or the drift of particles that can make the technique difficult to apply. A number of recent improvements proposed by the computer vision community [START_REF] Okuma | A boosted particle filter: Multitarget detection and tracking[END_REF][START_REF] Vermaak | Monte Carlo filtering for multi-target tracking and data association[END_REF][START_REF] Khan | MCMC-based particle filtering for tracking a variable number of interacting targets[END_REF] or specific to the biological field, like the introduction of multiple models of motion [START_REF] Smal | Particle filtering for multiple object tracking in dynamic fluorescence microscopy images: Application to microtubule growth analysis[END_REF] have made the method more suitable for biological tracking applications.

The Bayesian framework is still very attractive for particle tracking in biology, as exemplified by recent works [36,37]. In [36] Bayesian filtering was combined with a post processing technique to link incomplete tracks and obtain very long 2D trajectories of mRNA molecules in cells. In [37] the Bayesian tracking approach was used to study the in vivo colocalization of Oskar mRNA and trans-acting proteins.

While they have allowed significant improvements in the modeling of multiple particle tracking systems, the described Bayesian approaches do however suffer from a number of limitations in biological imaging. In particular, as they fully depend on the kinetic models of targets, which in many real situations are not discriminant enough or lack robustness, the association performances are often sub-optimal. For instance, in the frequent case of a high density of similar and closely spaced targets, subjected to Brownian diffusion, the association scores based on kinetic parameters are often not able to provide the true associations. Indeed, because of the randomness of the motion process, many hypotheses in the association problem have similar kinetic scores and it often happens that the true association hypothesis is not the one with the best kinetic score. In such cases, decisions based only on the kinetic information therefore lead to false associations that have a dramatic impact on the downstream track analysis.

One proposed solution to this problem has been to use the intensity of the detected objects in the Kalman filtering technique [START_REF] Smal | Particle filtering for multiple object tracking in dynamic fluorescence microscopy images: Application to microtubule growth analysis[END_REF]36,[START_REF] Genovesio | Multiple particle tracking in 3D+T microscopy: methods and application to the tracking of endocytosed quantum dots[END_REF]. The state space representation allows the extension of the measurement and state vector with a number of non kinetic features. The main difficulty of this kind of approach is to fuse this additional feature information with the kinetic information in a proper way. A probabilistic model of temporal evolution and a measurement model for non kinetic features are indeed required, and it is a major issue to design a balanced score function of the kinetic and the non kinetic information. Some of our recent work [START_REF] Chenouard | Feature-aided particle tracking[END_REF][START_REF] Chenouard | Improving 3d tracking in microscopy by joint estimation of kinetic and image models[END_REF] proposes a solution to these issues by integrating image and kinetic information in a statistical framework.

Joint image and kinetic based tracking

Tracking numerous and closely spaced targets is a very challenging task because of the limitations of the kinetic information and the interactions between targets. In order to solve this issue we have presented in [START_REF] Chenouard | Feature-aided particle tracking[END_REF] a unified Bayesian framework that fully exploits all the information available about the scene. We revisited the original definition of the scene likelihood, in order to define an enriched likelihood of association that combines in one function the target motion and the image formation process description, hereafter called appearance information.

The kinetic information and appearance information are split following the Bayes rule.

On one hand the kinetic part is computed thanks to a Kalman filtering-derived technique and using multiple models of motion. On the other hand the appearance likelihood relies on the description of a statistical model of image generation, which depends on the acquisition device. In practice this process is known accurately since the optical study of microscopy devices is a very active research field and accurate models exist [START_REF] Mahajan | Optical imaging and aberrations: Part I, Ray geometrical optics[END_REF][START_REF]Optical imaging and aberrations: Part II, Wave diffraction optics[END_REF]. We describe an image at a given time as the sum of the intensities originating from the targets, of a locally constant background and of a random acquisition noise. The noise is a random value that follows a probabilistic distribution that is generally well characterized, so its distribution will be used to give an accurate score to associations.

The difference between the observed image and the profiles and the background, which we call the residual, should follow the same distribution as the noise if the association is correct. Hence we use the statistics of the residual to compute the likelihood of an association given the image.

The appearance model is also used to split closely spaced targets that appear fused due to the resolution limit. The procedure is naturally integrated in the tracking procedure without any ad hoc decision step. Super resolution techniques [START_REF] Thomann | Automatic fluorescent tag detection in 3D with super-resolution: application to the analysis of chromosome movement[END_REF][START_REF] Thomann | Automatic fluorescent tag localization II: improvement in super-resolution by relative tracking[END_REF] also exploit the appearance of objects to address this issue, although the approach introduced in [START_REF] Chenouard | Feature-aided particle tracking[END_REF] performs significantly better, especially when the level of noise is high. This higher sensibility and robustness of the method comes from the exploitation of all the information available. Data such as dynamics of targets, statistics of targets appearance and disappearance, are used to decide if spots have to be separated. This method combined with the IMM filter was used to track in 3D hundreds of small fluorescent beads diffusing in water [START_REF] Chenouard | Improving 3d tracking in microscopy by joint estimation of kinetic and image models[END_REF]. Figure 1 depicts the trajectories extracted from this very dense environment where beads have a fast and directed motion due to a global flow. It also shows that the separation of aggregated particles is efficient.

Motion analysis and applications

Trajectories analysis

The final stage consists in the analysis of the trajectories resulting from cell or particle tracking. The exact type of analysis depends intrinsically on the biological question, but the study of the diffusion is most of the time a relevant issue. The Mean Square Displacement (MSD) [START_REF] Qian | Single particle tracking analysis of diffusion and flow in two-dimensional systems[END_REF][START_REF] Saxton | Single-particle tracking: Applications to membrane dynamics[END_REF] is a useful measure to characterize the diffusion of individual particles. By definition, the MSD is a function of the time lag during which displacements are measured and averaged. The shape of the MSD-time curve for a given trajectory is indicative of the mode of motion of the corresponding particle. We summarize the behavior of the MSD curve for different types of motion in Figure 2. In the case of normal diffusion by pure thermally driven Brownian motion, the MSD will increase linearly as a function of time. The diffusion constant determines the slope of the MSD line. In the case of flow or active transport the MSD increases more rapidly: quadratically with the time lag instead of linearly. The case of anomalous subdiffusion, characterized by a lagging MSD-time curve compared to normal diffusion, occurs if the motion is obstructed by obstacles. Confined motion, caused by corrals, tethering or other restrictions, manifests itself by a converging curve, where the limiting MSD value is proportional to the size of the region accessible for diffusion. Once the type of motion is identified, some subsequent parameters can be computed by regression. For instance the Brownian motion is fully characterized by its diffusion coefficient, and a directed motion is described by its direction and the velocity along it.

3D tracking of viruses

In [START_REF] Arhel | Quantitative four-dimensional tracking of cytoplasmic and nuclear HIV-1 complexes[END_REF], 3D tracking MSD analyses were performed in order to study virus-host cell interactions during infection of human cells by HIV-1 virus. The integrase of the virus was labeled with the bis-arsenical fluorescein derivative FlAsH which allowed imaging A recent study [START_REF] Gousset | Prions hijack tunnelling nanotubes for intercellular spread[END_REF] used latest advances in particle tracking methods [START_REF] Chenouard | Morphological source separation for particle tracking in complex biological environment[END_REF] to characterize the intercellular spreading of prions. Imaging techniques were used to show that TNTs allow transfer of exogenous and endogenous prions between infected and naive neuronal CAD cells. The CAD cells were labeled with LysoTracker red and imaged by spinningdisk confocal microscopy in 3D. Lyso-Tracker-positive vesicles were visualized moving towards naive cells inside nanotubes. The infectious vesicles were automatically detected, and their trajectories were reconstructed thanks to a Bayesian tracking method [START_REF] Chenouard | Morphological source separation for particle tracking in complex biological environment[END_REF], while moving inside TNTs and entering into the cells. A typical 3D prion vesicle trajectory is shown in Figure 3. The MSD analysis of the track revealed that the vesicle has a directed motion towards the recipient cell. The estimated velocity of the tracked vesicle (41.5 nm/s on average) is very similar to that calculated for a complex myosin VI and the GLUT1 transporter binding protein, indicating the possible involvement of actin

  mechanisms. For instance, viral infection processes include many steps involving dynamic interactions with different cellular structures such as the binding of the virus to specific receptors or attachment factors on the cell surface, and the transport to specific sites mediated by the cell cytoskeleton. The diversity of these interactions, and the possible existence of alternative infection routes, makes it necessary to study individual viral particles paths. Single particle tracking allows one to visualize and characterize the path of dynamic nano-metric objects, and provides a straightforward and accurate means to understand pathogenic mechanisms. Single-virus tracking has consequently emerged over the past few years as a powerful technique for the study of virus trafficking (see [3] and references therein for a comprehensive review) or prions spreading [4]. This review describes a number of particle tracking methods in time-lapse microscopy sequences, from the earlier and generally simpler techniques up to the most recent and more evolved tracking methods. It also provides recent examples from our own work that used particle tracking algorithms to investigate infectious diseases mechanisms. 2 Particle tracking in fluorescent microscopy 2.

  be tracked only if they are sufficiently distant and can be considered as spatially separated.

  intranuclear HIV-1 complexes. Tracking experiments of virus complexes resulted in the construction of 923 trajectories through time in 3D. The MSD analysis was used to determine the types of movement of HIV-1 complexes in the cytoplasm, during the association with the nuclear membranes and in the nucleus. For cytoplasmic HIV-1 complexes at early time points after infection, the shape of the MSD curve was a polynomial of order 2 which indicates a directed motion. Two types of such directed motions were distinguished based on their characteristics and velocities. First some fast and curvilinear motions (82 +/-7.4 nm/s) exhibited a saltatory property: both anterograde and retrograde movements were observed but with a global directionality towards the nucleus. This feature is compatible with microtubule-associated movement.The second class of motion was observed in the perinuclear area. It was still a directed motion but with lower speed (18+/-1.3 nm/s) on short distances, which is typical of actin-based transport. The implication of microtubules and actin in the motion of virus complexes was then confirmed by other standard biological techniques. For instance colocalization studies showed that HIV-1 complexes overlapped with both microtubules and actin filaments. Then pharmacological experiments were done: the dynein-based transport was uncoupled with a specific drug which resulted in the inhibition of entrance of the labeled virus into the nuclear compartment, while the inhibition of the actinmediated transport resulted on a drastic diminution of cytoplasmic transport of viruses.The dynamic study of association of complexes with nuclear membranes revealed significantly different characteristics compared to cytoplasmic transport. For instance, once the virus complexes were associated with the membrane their MSD trajectory curve was indicative of a restricted motion which implied confinement to a small volume or tethering.3.3 Prion tracking inside tunnelling nanotubes (TNTs)

mediated motors 4

 4 ConclusionRecent advances in optical systems and imaging protocols have allowed observing pathogens over time at the cellular and sub-cellular scale in multiple colors and dimensions in a fast and automated way. Systematic quantification of the recorded data sets requires fully automated image processing tools and algorithms specifically adapted to the problem at hand. Throughout this review we have browsed a number of methods developed to track sub-cellular objects in image sequences. The particle tracking problem includes many issues such as particles localization, motion estimation and measurements linking. The definition of the Bayesian framework has been an important achievement in the modeling of such a complex problem since it allows one to include in a uniform way prior knowledge about the scene and to define general objective functions. Enriching this framework with even more models dedicated to bioimaging and to solve problems in this formalism should attract much interest in the future. Nicolas, E., Chenouard, N., Olivo-Marin, J.-C., Guichet, A., A dual role for actin and microtubule cytoskeleton in the transport of golgi units from the nurse cells to the oocyte across ring canals. Mol.Biol. Cell 2009, 20, 556-568.

  Jaqaman, K., Loerke, D., Mettlen, M., Kuwata, H. et al., Robust single-particle tracking in live-cell time-lapse sequences. Nat. Methods 2008, 5, 695-702. [37] Mhlanga, M. M., Bratu, D., Genovesio, A., Rybarska, A. et al., In vivo colocalisation of oskar mRNA and trans-acting proteins revealed by quantitative imaging of the drosophila oocyte. PLoS ONE 2009. In press.

Figure 1 Fluorescent

 1 Figure 1

Figure 2 Mean

 2 Figure 2

Figure 3 Prion

 3 Figure 3

Wiley-VCH Biotechnology Journal


Dr Olaf Selchow
email: selchow@zeiss.de

Bernhard Goetze

Carl Zeiss

Productmanagement Microimaging Gmbh

Biosciences Goetze

Dr Bernhard Götze

Keywords: confocal, laser scanning microscope, variable secondary dichroic

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The new LSM 700 from Carl Zeiss.

The new LSM 700 confocal microscope is a member of the seventh generation of confocal microscopes from Carl Zeiss -a product family that is characterized by a wealth of genuinely innovative ideas and technologies. It is perfectly suited for a range of established advanced fluorescence microscopy applications like FRAP (Fluorescence Recovery after Photobleaching), FRET (Förster Resonance Energy Transfer) and Spectral Imaging. At the same time it also serves as a reliable and most sensitive workhorse for routine multicolor fluorescence and live-cell confocal imaging. The LSM 700 concept combines Carl Zeiss quality and exceptional ease of operation at an excellent price-performance ratio.

It can be equipped with up to 4 diode lasers (405 or 445, 488, 555, 639 nm) to cover the excitation spectra of the vast majority of commonly used fluorescence labels. Inside the scanning module the LSM 700 hardware is based on the same high end components as the Carl Zeiss Premium Class confocal microscope LSM 710. Examples are the patented PCT Laser coupling concept between laser fibers and scanning module, the beam combining optics as well as the solid state pinhole mechanics. Also scanners, PMT detectors, detection electronics and the Fixed Gate main beam splitter match the high end standards of Carl Zeiss premium systems. Together with a straightforward and impressively light-efficient beam path these components ensure the extraordinary detection sensitivity of the LSM 700. Even demanding applications with faintly stained fluorescent samples can be imaged with great image quality using only minimal excitation laser power (Figure 1).

The LSM 700 is available in a one or two channel version and the scanning module can be fitted to almost any Carl Zeiss high end scientific microscope. The system is operated by the next generation of ZEN, the user-friendly and popular Carl Zeiss LSM software platform. The ergonomically designed user interface as well as dedicated features like, e.g., the Smart Setup tool contributes to make the LSM 700 an extremely easy to use microscope system. Altogether, the LSM 700 is geared to cover all standard LSM applications with an outstanding performance. At the heart of the LSM 700, the new "Variable Secondary Dichroic" (VSD) splits the emission between two detectors (in the 2 channel configuration only). The VSD is a substrate coated with a special dichroic gradient layer. Embedded in a new detection beam path design the VSD is mounted on a motorized slider, and depending on its position, it can be a short pass/long reflectance dichroic at any wavelength between 420 and 630 nm (Figure 3). Thus, the user of the LSM 700 can freely choose the exact wavelength for the dichroic split. For example: With the VSD set to 527 nm, all the wavelengths shorter than 527nm will pass to channel 1, and all the wavelengths longer than 527 nm will be reflected to channel 2. The full integration of the VSD in the ZEN 2009 software allows switching the split from one scan to the next to ensure optimal fluorescence signal separation in the 2 detectors (Figure 4). Combining VSD position and emission filters, the user can freely adjust the active spectral detection window. The LSM 700 is thus prepared for any fluorescent dye combination and with the line-switching scanning mode 4 dyes can easily be imaged simultaneously at liveimaging frame rates. In addition to optimized channel separation, this unique detection beam path design also guarantees highest efficiency with regard to photon yield: No photon is ever lost since all the signal is directed either to detector 1 or 2.

Moved in defined steps between each scan, the VSD further permits the acquisition of spectrally resolved image series known as lambda stacks. These are prerequisite for getting optimal results in the separation of overlapping fluorescence signals by means of spectral unmixing. This technique allows the separation of even strongly overlapping fluorescent spectra like the signals from, e.g., GFP and YFP. Due to the unique detection beam path design and the properties of the VSD, each pair of images of a lambda stack represents the complete emission light. The spectral information is encoded in the relative changes of the image intensities with respect to the splitting wavelength. This property of the LSM 700 is the basis of a unique new way of lamda stack image acquisition and spectral unmixing.

The general theme of the development of the LSM 700 was to provide the user with a confocal microscope at no compromise in image quality and reliability: "Whatever this system does it does it extremely well". The technical concept to guarantee outstanding image quality, reliability and reproducibility of the acquired data includes, for example, an internal laser power calibration that guarantees constant light exposure during long time-lapse imaging or when re-using imaging parameters after months or years. Furthermore, a user-accessible software tool for automated calibration and alignment helps ensuring permanent perfect alignment and informs the user about the status of the system.

In summary, Carl Zeiss has set new standards by proving that uncompromised image quality and ease of use can come at an attractive price tag. The LSM 700 is a small, extremely sensitive and reliable versatile confocal microscope system which covers a vast range of applications like FRAP, FRET, spectral imaging, live-cell imaging and many more.

  Figure captions:

Figure 1 :Figure 4 :

 14 Figure 1:

Figure 5 :

 5 Figure 5: The LSM 700 can be fitted to almost any Carl Zeiss high end scientific microscope.

Wiley-VCH Biotechnology Journal