Francois Carlotti

Plankton

Keywords:

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

GLOSSARY

SPM : Structured population models IBM : Individual-based models ROMS: Regional Ocean Modeling System KtW model: Killing-the-Winner model Darwin model: high-resolution ocean and ecosystem model with many plankton species ECOPATH: free ecological/ecosystem modeling software (www.ecopath.org)

INTRODUCTION

The general purpose of population models of plankton species is to describe and eventually predict the changes in abundance, distribution, and production of targeted populations in their habitat conditions, i.e. under forcing of the abiotic environment, food conditions and predation. Computerbased approaches in plankton ecology were introduced during the 1970 with the application of population models to investigate large-scale population phenomena by the use of mathematical models.

Today, virtually every major scientific research project of population ecology has a modeling component. Population models are built for four main objectives, eventually for some practical purposes of management: (1) to estimate the survival of individuals and the persistence of populations in their habitat, and to look at the factors and processes which regulate their variability,

(2) to estimate the flow of energy and matter through a given population, and its role in trophic food webs;

(3) to study different aspects of behavioral ecology, and (4) to evaluate how genetic heterogeneity among populations impact the population structure of species and to identify transformations in biological traits that provide the adaptation to specific habitat condition.

The study of interactions between populations and successions of population are also topics related to population models. The field of biological modeling has diversified, and at present, complex mathematical approaches like neural networks, genetic algorithms and dynamical optimization start to be used thanks to the availability of supercomputers. However, the use of models in marine research should always be accompanied by extensive field data and laboratory experiments, either for initialization, verification or falsification, or continuous updating.

In recent years, plankton ecology has significantly grown from new data types with increased resolution of data acquisition with video, optical and acoustic plankton sensors on various platforms in relation with the physical environment, as well as from more comprehensive studies of organism physiology, behavior and life-history traits. New approaches to model marine population have 6 emerged (i.e. trait-based and size-based modeling) using more mechanistic descriptions of organism biology improving its modeled representation (Lindeman et al 2017).

APPROACH FOR MODELLING PLANKTON POPULATION

Population structure and units

A population is defined as a group of living organisms all of one species restricted in a given area and with limited exchanges of individuals from other populations. Thus, gene flows between quasi isolated populations are low in comparison with the significant gene flows within each population.

The first step in building a population model is to identify state variables (components of the population) and to describe the interactions between these state variables and external variables of the system and among the components themselves. The components of a population can be: 1) the entire population (1 component); 2) groups of individuals identified by a certain state: developmental stages, weight or size classes, age classes (fixed numbers of components), 3) all individuals or cohorts of individuals of similar history (varying numbers of components).

The usual unit in population dynamics models is the number of individuals per volume of water, but the population biomass can also be used (in g biomass or C or N). When all individuals or groups of individuals are represented, the individual weight can be considered as a state variable.

The forcing variables influencing the population dynamics are biological factors, mainly trophic nutriments and predators, and physical factors, mainly temperature, light, advection and diffusion.

Individual and demographic processes

Population models usually include four major processes: individual growth, development, reproduction and mortality. Growth is computed by the rate of individual body weight change.

Development is represented by the change of states (phases in phytoplankton and protozoan cell division, developmental stages in zooplankton) through which each individual progresses to reach maturity. Reproduction is represented by the production of new individuals. Mortality induces loss of individuals, and can be divided in two components: natural physiological mortality (due for instance to starvation) and predation. Combination of these four processes permits one to simulate (1) increase in terms of number of individuals in the population, (2) the body growth of these individuals, and (3) by combination of the two previous simulated values, the increase in total biomass of the population (which is usually termed "population growth").

Plankton Characteristics

The marine plankton community comprises a very large range of organisms from virus and prokaryotes to eukaryotes, uni-and multicellular, invertebrates and vertebrates, covering a wide size range and presenting a great variety of life cycles. Therefore the initial conceptual framework to build a population model structure requires identification of scales of the main biological characteristics of the studied species (Figure 1). Thus, essential information to build models of plankton organisms are: (1) the individual life duration (few hours for bacteria, one to few days for phytoplankton and unicellular animals, several weeks to years for zooplankton and ichthyoplankton organisms) and stage durations for those with complex life-cycles, (2) the range of change in size or weight between the beginning and the end of a life cycle (which is limited to a doubling in size in unicellular organisms, whereas the change in size and weight for metazoans are of several orders of magnitude and through different life stages), (3) the behavior and the swimming capacity on vertical and horizontal axis; and (4) the number of individuals produced by a mother individual (two individuals up to thousands of individuals). The observation time step has to be defined to adequately follow the time scale of the chosen variables, and thus should be lower than the duration of the shortest phases.

Plankton population models

Modeling plankton dynamics coupled to physical models and nutrient dynamics allows advancing our understanding of how pelagic marine ecosystems respond to climate variability. Therefore, the most modeled component in marine planktonic ecosystems is phytoplankton production. Most of the earliest phytoplankton models simulated the growth of phytoplankton as a whole, using only the process of photosynthesis. Nowadays, biogeochemical models consider different functional groups of phytoplanktonic populations, and a few models are dedicated to phytoplankton population growth dynamics at the species level including functional traits. Existing models of other unicellular plankton organisms (bacterioplankton, microzooplankton) usually treat them as a single unit, except for a few models simulating microbial cell cycles. On the contrary, mesozooplankton organisms, including the planktonic larval stages of benthic species (meroplankton) and fish with complex life cycles are extensively modeled at the population level.

Several types of modeling approaches can represent population demographics, including physiological traits and intra-population variability.

DYNAMICS OF SINGLE SPECIES

Populations models described by the total density

When a population is observed at timescales much larger than the individual life span, and on a large number of generations, models with one variable (the total number of individuals or the total biomass in that population) are the simplest. These models postulate that the rate of change of the population number, N, is proportional to N (eqn [1], where r is the difference between birth and death rates).

dN dt rN = [1]
The logistic equation ([2]) represents limitation due to the resources or space (see Figure 2).

dN dt rN N K = - ( ) 1 [2]
where K is the carrying capacity.

Population growth of bacteria, phytoplankton, or microzooplankton can be simulated adequately by the logistic equation. With addition of a time delay term into the logistic equation, oscillations of the population can be represented.

Figure 2 The growth of a plankton population with density regulation and its fitting by the logistic equation.

Population models of organisms with description of the life cycle

If observations of a population are made with a time step shorter than the life cycle duration (Figure 3), the population development pattern is a succession of periods with decreasing abundance of individuals due to mortality, and increasing abundance due to recruitment of new individuals in periods of reproduction (cell division or egg production). Recruitment is defined as the input flux of individuals in a given state (stage, size class, etc.).

Number of individuals Time

Production of new individuals

Loss by mortality To represent such patterns, it is necessary to identify different phases in the life cycle, based on age, size, developmental stages, and so on. Two types of models can be developed:

Life cycle

-Structured population models (SPM) which consider the flux of individuals through different classes.

-Individual-based models (IBM) which simulate birth, growth and development through stages and death of each individual.

These models are particularly used for representing the complex life cycles of zooplankton and ichthyoplankton, but have also been useful for studying the population growth of bacteria, phytoplankton, and microzooplankton, particularly division synchrony in controlled conditions Matrix models constitute a class of population models that incorporate some degree of individual variability. Matrix models are powerful tools for analyzing, for example, the impact of life history characteristics on population dynamics, the influence of current population state on its growth potential, and the sensitivity of the population dynamics to quantitative changes in vital rates. Matrix models are convenient for cases where there are discrete pulses of reproduction, but not for populations with continuous reproduction. They are not suitable for studying the dynamics of populations that live in fluctuating environments. Matrix population models [START_REF] Caswell | Matrix Population Models: Construction, Analysis, and Interpretation[END_REF]

 = + = m j 1 1 t 1, t j, N j F N [4]
The system of eqns [3] and [4] can be written in matrix form (eqn [5]).

( ) ( ) t N m N N N G m G G F m F F t N m N N N                                 - = +                           3 2 1 0 1 0 0 0 0 2 0 0 0 0 1 3 2 0 1 3 2 1 [5]
The Leslie matrix can easily be modified to deal with size classes, weight classes and developmental stages as the key individual characteristics of the population. Organisms grow through a given stage or size/weight class for a given duration.

There are several variations of matrix models, differing mainly in the expression of vital rates, which can vary with time depending on external factors (e.g., temperature, food concentration, competitors, predators) or internal (e.g., density-dependent) factors. For instance, the Lefkovitch matrix is a stage-structured or size-structured model of population growth, used when demography depends on the physiological stage. The transition from one stage to another is not a time function identical for all individuals (age). It is a function depending on size or weight allowing for the inclusion of intraspecific variation in growth rates (i.e. .due to different genotypes). There are many development using matrices to represent processes of multiple populations including genetic models, dynamics of metapopulations or multispecies models [START_REF] Newman | Modelling population dynamics: model formulation, fitting and assessment using state-space methods[END_REF] The earlier type of continuous-time structured model is usually referred to as the McKendrick-von Foerster equation, and uses the age distribution on a continuous-time basis in partial differential equations. This type of model has been developed to the extent that it can be used to describe population dynamics in a fluctuating environment. In addition, it also applies to situations in which more than one physiological trait of the individuals (e.g., age, size, weight and energy reserves) have strong influences on individual reproduction and mortality. The movement of individuals through the different structural classes is followed over time. Age and weight are continuous variables whereas stage is a discrete variable.

The general equation is eqn [6], where n is abundance of individuals of age a and mass m at time t.

m) a, n(t, m) a, (t, m m) a, m)n(t, a, (t, a m) a, n(t, t m) a, n(t,        - = + + g [6]
The von Foerster equation describes population processes in terms of continuous age and time (agestructured models) according to eqn [7].

     n(a, t) t n(a, t) a (a, t) n(a, t) + = - [7]
The equation has both an initial age structure  at t=0: (eqn [8]) and a boundary condition of egg production at a = 0 (eqn [9]).

n(a,0)

0 (a) =  [8] n(0, t) F(a,S 0 ) n(a, t) da R =   [9]
F is a fecundity function that depends on age (a) and the sex ratio of the population SR. These kinds of equations are mathematically and computationally difficult to analyze, especially if the environment is not constant.

The same type of equation as [7] can be used where the age is replaced by the weight (weightstructured models, substituting w to a in equation7). The weight of the individual w and the growth g are influenced by the temperature T, the food P and by the weight itself through allometric metabolic relationships, and F is the fecundity function which depends on weight (w) and the sex ratio of the population SR. The numerical realization of this equation requires a representation of the continuous distribution n(w,t) by a set of discrete values ni(t) that are spaced along the weight axis at intervals wi = wi+1 -wi.

Stage-structured population models

Plankton populations often have continuous recruitment and are followed in the Reld by observing stage abundances over time. A large number of zooplankton population models deal with population structures in term of developmental stage, using ordinary differential equations (ODEs).

16

A single ODE can be used to model each development stage or group of stages: for instance, a copepod population can be subdivided into four groups: eggs, nauplii, copepodites, and adults. The equation system is eqns [10]- [13], where R is recruitment,  is the transfer rate to next stage, and  is the mortality rate. The system of ODEs is solved by Euler or Runge-Kutta numerical integration methods, usually with a short time step (approximately 1hour).

eggs dN dt R N N 1 1 1 1 1 = - -   [10] nauplii dN dt N N N 2 1 1 2 2 2 2 = - -    [11] copepodids dN dt N N N 3 2 2 3 3 3 3 = - -    [12] adults dN dt N N 4 3 3 4 4 = -   [13]
In the model presented in Mar. Ecol. Prog. Ser.,56,3,[225][226][227][228][229][230][231][232][233][234][235][236][237][238][239][240][241][242].

Individual-based models of a population

Because transport of planktonic organisms is largely controlled by oceanic currents, numerous IBM of planktonic populations, or of planktonic larval forms of nektonic and benthic populations, have been developed and coupled to hydrodynamic models of varying degrees of three-dimensional and temporal complexity. IBMs have been applied for phytoplankton and zooplankton species, meroplanktonic larvae and early life-history of fish populations. Individual-based modeling approach is well established for metazoans where organisms are more functionally complex, but are increasingly being applied to unicellular organisms when growth kinetics are formulated as a function of internal properties. Object-oriented programming (OOP) and cellular automata are techniques that have been applied to IBMs recently.

IBMs do not make the same averaging assumption as SPMs, which can be important when model equations are nonlinear and heterogeneity exists, and they are more suitable for simulating complex The model starts with an initial population and the basic environment, then monitors the changes of each individual. At any time t, the i-state of individual j changes as eqn [14].

Xi,j(t) = Xi,j (t-dt) + f(X1,j (t-dt),...., Xi,j (t-dt) [14] Xi,j (t) is the value of the i-state of individual j, and f is the process modifying Xi,j, as a function of the values of different i-states of the organism, external parameters such as the temperature.

When the fate of all individuals during the time-step dt has been calculated, the changes to the environment under the effects of individuals can be updated. Any stochastic process can be added to eqn [14].

This type of model can add a lot of detail in the representation of physiological functions. Individual growth can be calculated as assimilation less metabolic loss, and the inter-individual variation in physiology can be represented by adding stochastic processes or parameters describing the characteristics of each individual (growth and development parameters, mortality coefficient, and parameters connected with reproduction). The end results are unique life histories, which when considered as a whole give rise to growth/size distributions that provide a measure of the state of the population Individual traits and population structure and dynamics emerge from the properties of metabolic organization and the indirect interactions of individuals via competition for food.

Calibration of parameters

Parameterization of a model can range from very simplistic to extremely complex depending upon the amount of information known about the population under consideration. Bioenergetic processes

Spatial plankton dynamics with advection-diffusion-reaction equations

Eulerian or grid-based models are used to simulate the distribution of variables that behave like tracers (e.g. plankton concentrations). The state of a variable in fixed space volumes (i.e. boxes or grid cells) is represented by a single average value and the dynamics in the field domain of the grid model depends on the biological changes within each box and exchanges between boxes. Equation

[15] is the general physical-biological model equation used to describe the interaction between physical mixing and biology.

( ) ( )

   C t C K C bio ical terms a +   -   =" log " [15] C(x,y,z,t,)
      z t) z, (x, C t) z, (x, w bi dynamics population z t) z, (x, C t) K(z, z x t) z, (x, C t) K(x, x z t) z, (x, C t) z, (x, w a x t) z, (x, C t) z, (x, u a t t) z, (x, C i i i i i i                 + =       -       - - - [16]
The population dynamics model was presented in equations [10] to [13]. The upwelling zone extended 50 km from the coast down to a depth of 50 m, and was divided into a grid with spacing 2.5 m in depth and 1 km in the horizontal. The author used a finite difference scheme with a time step of 1 hour, which fell within the bounds for computational stability.

Coupling IBMs and spatially explicit models

This approach uses simulated currents from sophisticated 3D hydrodynamic models driving Lagrangian models of particle trajectories to examine dispersion processes. IBMs of planktonic species are associated with sparsely distributed particles, each representing an individual planktonic organism or a group of identical individuals, which are forced with output from an Eulerian biophysical model (i.e. one-way coupling) of varying degrees of 3D and temporal complexity. The approach is relatively straightforward and is a first step in formulating spatially-explicit individual based models. Given a "properly resolved" flow field, particle trajectories are computed (generally with standard Runge-Kutta integration methods of the velocity field). Specifically, hydrodynamic models provide the velocity vector v = (u,v,w) as a function of location x = (x,y,z) and time t and the particle trajectories are obtained from the integration of eqn [17].

d x / dt = v(x,y,z,t) [17]
The simplest model of dispersion is a random walk model in which individuals move along a line from the same starting position. These trajectories could be modified by turbulent dispersion as described in the section below. Once the larval/particle position is known, additional local physical variables can be estimated along the particle's path: e.g., temperature, turbulence, light, etc., and input to the IBM. The physical quantities are then included in biological (physiological or behavioural) formulations of IBMs (see below).

By coupling an IBM population of key species with its environmental forcing (food, temperature)

provided by a coupled physical-biogeochemical model, the simulations allow to evaluate the influence of environmental variability of the bio-physical environment (flow, temperature, turbulence, light, and trophic resources) on dispersal, growth and mortality of planktonic life stages and to explain how patchy habitats may influence the distribution of individuals within a population and recruitment in a given area. Simulations considering trajectories of plankton as passive particles are usually a first step before considering any active swimming capability of planktonic organisms either using a prescribed vertical migration strategy, or integrating a vertical motion (sinking or swimming) which could depend on animal stage development and/or its physiological state. The simulated plankton distributions from these models tend to compare better with observed distributions than models that use passive particles. Sensitivity studies show that behaviour is an important factor in determining larval transport and/or retention.

The coupling of IBMs of zooplankton and fish populations and 3D circulation models is still a recent field of study (Figure7). Available on-line user-friendly design of lagrangian tools such as Ichthyop Generally, models that describe the spatial heterogeneity of the habitat have been designed to answer questions about the spatial and temporal distribution of a population, particularly about connectivity of metapopulations between habitats (e.g. Cowen et al 2009). They allow us to explore the potential effects of habitat alteration on these populations.

Modeling behavioral mechanisms, aggregation and schooling, patches

Different types of models have been built, some of them focusing on the structure and shape of aggregations depending on internal and external physical forces, others dealing with the benefits for individuals of living in groups with regard to feeding (foraging models) and to predation. The Lagrangian approach can take into account the behavior of individual organisms, and the effects of the physical environment upon them. Although Eulerian approaches are mathematically tractable, the methods do not explicitly address the density dependence of aggregating individual behavior within a patch.

Dynamic optimization allows descriptions of the internal state of individuals, which may lead to both variable and fluctuating motivations among individuals over short time periods

INTERACTIONS BETWEEN POPULATIONS

Models with plankton populations in interactions

Interactions between populations by competition for shared resources and by predation are among the most important factors driving plankton community dynamics. Interaction between marine species have been rarely treated by stage-resolved population models, although structured population models as well as IBM models can represent interactions between species like predation, parasitism or even cannibalism.

Simple models of two species interactions take the form of eqns [18] and [19] (The so-called "Lotka-Volterra equations" or predator-prey equations):

d N dt r N k N N 1 1 1 1 1 2 = - [18] d N dt r N k N N 2 2 2 2 1 2 = - [19]
Different extensions of the second term on the right part of the equation can be used. Such equations have been used in many situations of predator-prey interaction (e.g. bacteria-phage, carnivorous zooplankton -prey species) or between species of the same trophic level in competition for resources

As an example, Gaedke and Ebenhöh presented in 1991 (reference in [START_REF] Carlotti | Modelling zooplankton dynamics[END_REF] a study on the interaction between two estuarine species of copepods Acartia tonsa and Eurytemora affinis.

They first used a simple model based on eqns [18] and [19] representing two bulk populations in interaction, and secund two-stage-structured population models with similar eqns to [10]-[13])

including crossed interactions between stages of the two species. Only the second model resulted in stable coexistence between the two species due to the predation on larval stages of both species by Acartia tonsa.

The original "kill the winner" (KtW) model (see references in Lindeman et al. 2017) explaining the continual succession of different planktonic populations in competition for resources is formulated as deterministic Lotka-Volterra-type equations for the species biomass concentrations. Nowadays, KtW models have been applied to planktonic communities consisting of viruses; prokaryotes, phytoplankton and zooplankton protozoans and metazoans).

Food webs with population models

In the last decades, there has been an increasing use of ecosystem models including representation of multiple phytoplankton and zooplankton populations in interaction (Plankton Functional Type models) and coupled to physical models to study the effect of biodiversity on ecosystem functioning (e.g. the Darwin model, [START_REF] Follows | Emergent biogeography of microbial communities in a model ocean[END_REF]. The different species represent differing degrees of functional diversity defined by species' functional traits. Phytoplankton species compete for nutrient resources and represent resources for zooplankton species which additionally exert cannibalism, and so on for zooplankton predators. The planktonic community structure emerges from a range of possibilities determined by processes of dispersal, competition for resources, predation and physiological characteristics.

Although the amount of energy transferred through the pelagic food chain from phytoplankton to higher trophic levels can be represented by food web models taking into account bulk species biomass (Ecopath), size-based food-web approach is preferred by assuming species-and size-based interactions which result from the dominance of small primary producers, size-based predation and ontogenetic increases in the trophic levels of many consumers. The output of these models is that the community structure is an emergent property of the ecosystem functioning itself.

________________________________________________________________________________

Figure 1

 1 Figure 1 Schematic representation of key differences in biological rates, species numbers, and size/weight growth during the life cycle, and swimming behaviour between the different trophic levels in marine ecosystems. Modified from Carlotti F and Poggiale J.C. 2010. Towards methodological approaches to implement the zooplankton component in "end to end" food-web models. Progress in Oceanography 84: 20-38.

Figure 3

 3 Figure 3 Total population abundance in controlled conditions. The population is initiated with newborn individuals, and decreases, first owing to mortality, until the maturation period, when a new generation is produced. At the end of the recruitment period of new individuals of second generation, the total abundance decreases, again owing to mortality. Owing to individual variability in development, loss of synchronism in population induces a broader recruitment period for the third generation. In non-limiting food conditions, recruitment of new individuals is continued after few generations. The three generations correspond to the exponential phase presented in Figure 2 (i.e., without density regulation).

(

  chemostat). The major distinction between physiologically structured population models and individual-based models in a stricter sense is that individual-based models track the fate of all individuals separately over time including their states (e.g. body weight, reserve pool or cell quota) and behavior (e.g. swimming, predation), while physiologically structured population models follow population-level properties such as the density or biomass of individuals in specific classes (age or size classes, stages) using differential equations. For both model types, the physiological submodel can be represented with more or less details up to a complete budget model including a representation of structural growth and reserve pool, such as DEB models (Dynamic Energy Budget theory, Kooijman 2010)Structured population modelsA population can be structured with respect to age (age-structured population models), stage (stagestructured population models), or size or weight (size or weight-structured population models). Two types of equations systems are usually used: matrix models, which are discrete-time difference equation models, and continuous-time structured population models.

Figure 4 ,

 4 the transfer rate of animals from stage to stage and the mortality at each stage are expressed as simple linear functions, which induce a rapidly stable stage distribution. To represent the delay of growth within a stage, more refined models consider ageclasses within each stage, or systems of delay differential equations. They have a high degree of similarity with observed cohort development in mesocosms or closed areas (Figure5).

Figure 4

 4 Figure 4 Schematic representation of the population dynamics mathematically represented by eqn [10]-[13]. i : transfer rate of stage i to stage i+1; i: mortality rate in stage i ; R: recruitment: number of eggs produced by females per day.

Figure 5

 5 Figure 5 Simulation of the cohort development of the copepod Euterpina acutifrons in mesocosms with a structured stage and age-within-stage model. (A) Total population during development, with variable temperature and constant food supply (points = experimental data; line = simulation). The initial density decreases owing to mortality and then increases to new born individuals (as the first part in Figure 3). (B) Naupliar stages (N1 to N6), copepodite stages (C1 to C5), and adults during development (points = experimental data; line = simulation). The simulation starts with similar N1 of same age. Reproduced from Carlotti F and Sciandra A, 1989. Population dynamics model of Euterpina acutifrons (Copepoda: Harpacticoida) coupling individual growth and larval development.

  and adaptive behavior. IBMs treat populations as collections of individuals, with explicit rules governing individual biology and interactions with the environment. Each biological component can change as a function of the others. Each individual is represented by a set of variables that store its istate (e.g. age, size, weight, nutrient or reserve pool, etc.). These variables may be grouped together in some data structure that represents a single individual, or they may be collected into arrays (an array of all the ages of the individuals, an array of all the sizes of the individuals, etc.), in which case an individual is an index number in the set of arrays. The i-state of an individual changes as a function of the current i-state, the interactions with other individuals, and the state of the local environment. The local environment can include prey and predator organisms that do not warrant explicit representation as individuals in the model. Population-level phenomena (e.g., temporal or spatial dynamics) or vital rates can then be inferred directly from the contributions of individuals in the ensemble.

(

  ingestion, egestion, excretion, respiration and egg production) are usually modeled from experimental results, whereas biometrics (size, weight, …) and demographic (development rate, mortality rate, ) parameters are estimated by combining data from life table collected in the field or from laboratory studies.Calibration task involves systematic adjustment of model parameter estimates so that model outputs more accurately reflect observed data.. To solve for the unknown parameters, new techniques have been developed such as inverse methods and data assimilation by fitting simulations to data. The use of allometric relationships of physiological rates to weight or applying conceptual frame such as the DEB theory where underlying physiological processes based on first principles (e.g. massenergy conservation laws) allows reducing the number of parametersSPATIAL DISTRIBUTION OF SINGLE PLANKTON POPULATIONSAn important development in plankton population modeling is to make full use of the increased power of computers to simulate the dynamics of plankton (communities or populations) in sitespecific situations by coupling biological and transport models, giving high degrees of realism for interpreting plankton population growth, transport, spatial distribution, dispersion, and patchiness.These models fall into two basic categories to represent spatial plankton dynamics (1) Eulerian or grid-based models which follow the evolution of some quantity (abundance or biomass of structured population models) at discrete, fixed physical locations using advection-diffusion-reaction equations, and (2) Lagrangian biological models which follow individual organisms or groups of organisms moving through space as they are advected by the ambient currents with the biology of each organism described by an IBM module. Continuous development in physical -biological coupled modeling technique for both types of models have opened opportunities to study the variability of the seasonal population cycles under different scenarios of environmental forcing for different regions of the ocean[START_REF] Ji | Marine plankton phenology and life history in a changing climate: current research and future directions[END_REF]

  is the concentration of the biological variable at position x,y,z at time t (expressed as numbers of organisms or biomass of organisms per unit volume) which could be either a planktonic functional group (e.g. zooplankton), or a species or a size class (in which case the number of equations would equal the number of size classes). . a(u a,v a,w a) represents the advective fluid velocities in x,y,z directions. Kx, Ky, Kz are diffusivities in x,y,z directions.  = (/x, /y, /z) is the Laplacian operator.On the left hand side of eqn[15], the first term is the local change of C, the second term is advection caused by water currents, and the third term is the diffusion or redistribution term. The right hand side of eqn[15] has the biological terms that represent the sources and sinks of the biological variable at position x,y,z as a function of time. The biological terms may include a velocity component (swimming of organisms, migrations, sinking, ...), and the complexity of the biological representation can vary from the dispersion of one (the concentration of a cohort) to detailed population dynamics.Physical-biological models of various levels of sophistication have been developed recently for different regions of the ocean. Biological models can be configured as compartmental ecosystem models in an upper ocean mixed-layer, where phyto-, microzoo-and meso-zooplankton are represented by one variable. In extended case, the model takes into account several size classes of phyto-, microzoo-and zooplankton. Recent model developments consider concurrent species with variable functional traits (e.g. Darwin model,[START_REF] Follows | Emergent biogeography of microbial communities in a model ocean[END_REF]). Studies of plankton population distribution in heterogeneous marine landscape where plankton may be aggregated in dynamical mesoscale structures (e.g. upwelling-and downwelling-regions, Langmuir circulations, eddies) can be undertaken with population described by equations of the McKendrick-von Foester type coupled with 2D or 3D fine-grid hydrodynamical models(Runge et al., 2005).A classic example by Wroblewski in 1982 is the stage-structured population model of Calanus marshallae, a copepod species, embedded in a circulation system simulating the upwelling off the Oregon coast (Figure6). Simulations of the dynamics focused on the interaction between diel vertical migration and offshore surface transport. The zonal distribution of the life stage categories Ci of C. marshallae over the Oregon continental shelf was modeled by the two-dimensional (x,z,t) equation[16], where the vertical swimming speed of the ith stage wbi follows a sinusoidal function of time: wbi = wsi sin(2t), with wsi the maximum migration speed of the stage.

Figure 6

 6 Figure 6 Model of Calanus marshallae in the Oregon upwelling zone. The figure shows the simulated zonal distribution of (A) eggs, (B) nauplii, and (C) early copepodites at noon on 15 August. Concentrations of each stage are expressed as a fraction of the total population (all stages m -3 ). Reproduced from Wroblewski JS, 1982. Interaction of currents and vertical migration in maintaining Calanus marshallae in the Oregon upwelling zone-a simulation. Deep Sea Research 29: 665-686.

(

  http://www.ichthyop.org;Lett et al., 2008) allows to simulate zoo-and ichthyoplankton dynamics in any oceanic region by incorporating the most important processes involved in plankton or fish early life, and to force IBM module by input time series of velocity, temperature and salinity fields archived from ROMS.

Figure 7 :

 7 Figure 7: Simulated IBM of fish eggs and larvae dispersal pattern displayed on a longitude latitude plot with the ichthyop tool (http://www.ichthyop.org). Reproduced from Lett et al. (2008) A Lagrangian tool for modelling ichthyoplankton dynamics.. Environmental Modelling and Software, Elsevier, 23 (9), 1210-1214.

  See also: Biogeochemical Data Assimilation (xxx). Carbon Cycle (xxx). Fish Larvae (2xxx 6). FishPredation and Mortality (xxx). Fish Vertical Migration (xxx). Gelatinous Zooplankton (xxx). Krill (xxx). Lagrangian Biological Models (xxx). Large Marine Ecosystems (xxx). Marine Mesocosms(xxx). Microbial Loops (xxx). Models of Small-scale Patchiness (xxx). Network Analysis of Food Webs (xxx). Nitrogen Cycle (xxx). Ocean Gyre Ecosystems (xxx). Patch Dynamics (xxx). Phosphorous Cycle (xxx). Polar Ecosystems (xxx). Small-scale Physical Processes and Plankton Biology (xxx). Upwelling Ecosystems (xxx). Zooplankton (xxx).


Étienne Polack
email: etienne.polack@math.cnrs.fr

Yvon Maday
email: yvon.maday@sorbonne-universite.fr

Andreas Savin
email: andreas.savin@lct.jussieu.fr

A contribution to the

The Kohn-Sham method uses a single model system, and corrects it by a density functional the exact user friendly expression of which is not known and is replaced by an approximated, usable, model. We propose to use instead more than one model system, and use a greedy extrapolation method to correct the results of the model systems. Evidently, there is a higher price to pay for it. However, there are also gains: within the same paradigm, e.g., excited states and physical properties can be obtained. * Note however that this is at the prize of working in R 3N instead of R 3 , and thus requiring accurate many-body, e.g., configuration interaction calculations.

Introduction 1.Motivation

Density functional theory (DFT) has a weak point: its approximations (DFAs). First, the Hohenberg-Kohn theorem tells us that there is a density functional for electronic systems, F [ρ], that is universal (that is, indepen-1 dent of the potential of the nuclei), but does not give us a hint on how systematic approximations can be constructed. In practice, models are produced to be fast in computations, typically by transferring properties from other systems, like the uniform electron gas. Second, the most successful approximations are using the Kohn-Sham method (introducing a fermionic wave function) that decomposes F [ρ] into the kinetic energy, the Hartree energy and an exchange-correlation energy contribution although the question of how and what part of F [ρ] should be approximated is, in principle, open.

In the present contribution we totally change the paradigm in the following way still led by the issue of universality. Let us start with a physical consideration. When electrons are close, the Coulomb repulsion is so strong that some of its features dominate over the effect of the external potential. This is also reflected mathematically in the short-range behavior of the wave function, as present in the Kato cusp condition [START_REF] Kato | On the eigenfunctions of many-particle systems in quantum mechanics[END_REF][START_REF] Fournais | Analytic Structure of Many-Body Coulombic Wave Functions[END_REF][START_REF] Yserentant | Regularity and Approximability of Electronic Wave Functions[END_REF][START_REF] Flad | Explicit Green operators for quantum mechanical Hamiltonians. II. Edge-type singularities of the helium atom[END_REF], and in higher-order terms [START_REF] Rassolov | Behavior of electronic wave functions near cusps[END_REF][START_REF] Kurokawa | General coalescence conditions for the exact wave functions: Higher-order relations for coulombic and non-coulombic systems[END_REF]. We further note that approximating numerically the short-range part of the wave function needs special care, due to the singularity of the Coulomb interaction when the electrons are close.

The considerations above and the independence of the interaction between electrons from that between them and the external potential provides a basis for constructing approximations. Thus, we propose to solve accurately a Schrödinger equation with a Hamiltonian that is modified to eliminate the short-range part of the interaction between the electrons which is one of the difficult parts in the numerical simulations. The way to do it is not unique, and we try to turn this to our advantage: we use several models, and from them we try to extrapolate to the physical system [START_REF] Savin | Correcting model energies by numerically integrating along an adiabatic connection and a link to density functional approximations[END_REF]. In other words, we follow an "adiabatic connection" (see [START_REF] Harris | The surface energy of a bounded electron gassolid[END_REF]), without ever constructing a density functional. This new paradigm thus explores the possibility to replace the use of DFAs by mathematically controlled approximations: we make density functional theory "without density functionals."

Our approach has introduced an additional difficulty nonexistent in the Kohn-Sham method: the long-range part of the interaction has to be treated accurately, and not only its electrostatic component. One may ask whether this additional effort is justified, and whether one gains anything with respect to a calculation where the physical (Coulomb) interaction is used. For a single calculation, the gain is due to the lack of singularity in the interaction expressed by a weak interaction potential allowing for simplified treatments, such as perturbation theory. However, as the extrapolation to the physical system needs more than one point, it is essential that the number of points stays very small, and the interaction weak.

Objective and structure of the paper

We first choose, in Sec.2.1, a family of model (parameter dependent) Hamiltonians that are more flexible than using only the Kohn-Sham (noninteracting) Hamiltonian. * This is followed by a description of how universality is introduced, namely by analyzing how a nonsingular interaction approaches the Coulomb one, and not by transfer from other systems, as usually done in DFAs. The physical system of interest is one among the parameter dependent models corresponding to some precise value of the parameter; in Sec. 2.2 its solution is extrapolated from the solutions to the models for other values of the parameter, expected that these solutions are more simple to be approximated. This extrapolation is efficiently handled in the general framework of the model reduction methods and more precisely referring to a variation of the Empirical Interpolation Method [START_REF] Barrault | An 'empirical interpolation' method: Application to efficient reduced-basis discretization of partial differential equations[END_REF].

We believe that such an approach can not only discuss what DFAs are really doing, but can evolve to being used in applications. Some argument supporting this statement is given. However, in this paper numerical examples (gathered in Sec. 3) are only presented for two-electron systems that are numerically (and sometimes even analytically) easily accessible: the harmonium, the hydrogen anion, H -, and the hydrogen molecule, H 2 in the ground state, at the equilibrium distance.

As we do not use the Hohenberg-Kohn theorem, the technique can be applied without modification also to excited states. We provide in Sec. 3.5, as an example, the first excited state of the same symmetry as the ground state.

Some conclusions and perspectives are presented in Sec. 4. Finally, in order to facilitate reading the manuscript, various details are given in Appendices A-E that follow Sec. 4.

Approach

The model Schrödinger equation

We study a family of Schrödinger equations,

H(µ)Ψ(µ) = E(µ)Ψ(µ), (2.1)
where µ is some nonnegative parameter. More precisely, in this paper, we use

H(µ) = T + V + W (µ), (2.2)
where T is the operator for the kinetic energy, V is the external potential (in particular that of the interaction between nuclei and electrons) and W (µ) represents the interaction between electrons. Although not required by the general theory, in this paper we introduce the dependence on µ only by modifying the interaction between electrons,

W (µ) = i<j w(r ij , µ), (2.3)
choosing

w(r ij , µ) = erf(µr ij ) r ij (2.4)
where r ij = |r i -r j | is the distance between electron i (at position r i ) and electron j (at position r j ). Finally, the external potential V is written like

V = N i=1 v(r i ). (2.5)
where v is the local one particle operator. Note that the N -particle operators are denoted by upper case letters, while the one-particle operators are denoted by lower case letters. Note also that for µ = 0 we have a trivial noninteracting system, while for µ = ∞ we recover the Coulomb system. The operator w is long-ranged: as µ increases, the Coulomb interaction 1/r 12 starts being recovered from large distances. The first reason for this choice is that, as mentioned above, we expect a universal character for short range (this is related to the difficulty of common DFAs to correctly describe long-range contributions, cf. Appendix A). The second reason is that the solution of Eq. (2.1) is converging more rapidly with (conventional) basis set size when the interaction has no singularity at r 12 = 0.

In principle, introducing a dependence of the one-particle operators (T and V ) on µ makes the formulas a bit more clumsy, but does not introduce important difficulties in its application. Using such a dependence might improve the results, but it is not discussed in this contribution. In the following, in order to simplify notation, we drop the argument µ, when µ = ∞, e.g., E = E(µ = ∞).

2.2.

The correction to the model

Using a basis set

Of course, solving the Schrödinger equation for the model, Eq. (2.1) with finite µs, does not provide the desired solution, i.e., the one that is obtained for µ = ∞. We thus need to estimate the difference in eigenvalues:

Ē(µ) = E -E(µ).
(2.6)

Since Ē(µ) tends to zero at infinity, the idea is to first expand this difference Ē(µ) in a basis (of functions that tend to zero at infinity), retaining M terms,

Ē(µ) ≈ ĒM (µ) = M j=1 c j χ j (µ), (2.7)
leading to

E(µ) ≈ E - M j=1 c j χ j (µ),
or, more precisely since E is not known, we replace it by an approximation denoted as E M ,

E(µ) ≈ E M - M j=1 c j χ j (µ). (2.8)
The idea then proceeds by determining the unknown E M values and the coefficients c i from M + 1 values of E(µ m ), for m = 0, . . . , M for an appropriate choice of the parameter values µ m . Finally, taking into account that the functions χ j tend to zero at infinity, the proposed approximation for E is E M . Of course, this extrapolation approach often fails if care is not enough taken in the choice of the functions χ j , 1 ≤ j ≤ M , and the values µ m , for m = 0, . . . , M . First, one has to decide about their form. Second, one has to find a way to keep M as small as possible to reduce computational cost while preserving a good accuracy.

Approaching the Coulomb interaction

As recalled above, we derive from the leading term of the Coulomb interaction between the electrons that, to leading order, the solutions of differential equations are determined at short range by the singularities. The interaction w in Eq. (2.4) has no singularity at r 12 = 0, for any finite µ. However, as the parameter µ increases, w(•, µ) approaches the singular Coulomb potential.

In order to see how this limit is approached, let us perturb the exact solution. To first order, the perturbation correction to the energy is given by

Ē(µ) = Ψ W -W (µ) Ψ , for µ → ∞.
(2.9)

By changing the integration variables r i to µr i we see that

Ē(µ) ∝ µ -2 as µ → ∞, (2.10)
providing a leading behavior that we want the basis functions χ i to reflect. It is possible to continue this analysis for higher order terms. In fact, the next term (in µ -3 ) has a coefficient proportional to that of µ -2 , the proportionality coefficient being determined by the nature of the Coulomb singularity [START_REF] Goll | A short-range gradient-corrected spin density functional in combination with long-range coupled-cluster methods: Application to alkali-metal raregas dimers[END_REF].

Choice of the basis functions

In the main part of this contribution we use a simple ansatz,

χj (µ) = 1 -j µ(1 + j 2 µ 2 ) -1/2 , j = 1, . . . , M, (2.11)
that respects indeed the condition of Eq. (2.10). The motivation for this specific choice, that is arbitrary to a certain degree, as well as some results obtained with other choices of basis functions, is given in Appendix B.

The first functions of this basis set are presented in Fig. 1, together with an example of a function it has to approximate. It illustrates that the function we want to describe is between basis function χ2 -for small µand basis function χ3 -for large µ. However, a simple linear combination between these (only) two surrounding basis functions from the family in Eq. (2.11) does not improve much the accuracy, but of course, more (and more appropriate) functions in the family can (and should) be called.

Reducing the basis set

Using a large set of χ j (a large M ) can rapidly become computationally prohibitive (because it requires a large number of evaluations of E(µ m ), for m = 0, . . . , M ) and numerically unstable (because it is classically much more difficult to stabilize extrapolation than interpolation). In order to reduce their number and increase the stability of the extrapolation, we use a greedy (iterative) method, as in the Empirical Interpolation Method (EIM) leading to proper choices of µ m , for m = 0, . . . , M known as "magic points."

In the Kth iteration of EIM, one starts from a set of K-1 basis functions (for us, χj ) and K -1 points (for us, μj belonging to some (discretized) interval, say, close to zero, to benefit at most of the regularization of the erf function). One then chooses the Kth function χK-1 (among the remaining M -K basis functions) as being the one that is most poorly approximated by the current interpolation (based on the K -1 basis functions and the K -1 points) in a sense dedicated to the final goal we want to achieve (that can be uniform error, error on some part of the domain, or even at some value) and the Kth point μK-1 that, in the admissible set, brings the more information. In this contribution, as we are only interested at extrapolating the value of µ at infinity, so we chose the error as the absolute value of difference between the Kth basis function and its interpolant at infinity as the final goal we want to achieve.

Note that the procedure selecting the next point and function does not make any use of the function to be approximated (here Ē). It is thus a cheap step compared with the calculation of E(μ m ) on the system of interest.

To improve the results for the extrapolation, we have modified the EIM algorithm into what we call the Forward Looking EIM (FLEIM). While EIM tries to get the maximal improvement through a sequential choice of, first the new basis function, then the new point of interpolation, FLEIM tries to get the best pair for improvement in the selected goal. The method is explained in more detail in Appendix C. In what follows, we present the results of FLEIM as they are better and more stable than those of EIM, as is illustrated in App D.1.

Computing other physical properties

FLEIM can be used to approximate other physical properties, i.e., correct expectation values of operators A = H obtained with the model wave functions, Ψ(µ),

A(µ) = Ψ(µ)|A|Ψ(µ) .
(2.12)

This can be seen immediately by noting that the derivation in Sec. 2.2.1 is not specific for correcting E(µ), but can also be applied to A(µ).

For the choice of the basis functions, we point out that properties are obtained by perturbing the Hamiltonian with the appropriate operator, say, A, H → H(λ) = H + λA.

(2.13)

The expectation value of A can be obtained as the derivative of E(λ) w.r.t. λ, at λ = 0. Of course, this procedure can be applied to model Hamiltonians, yielding E(λ, µ) and

Ψ(µ)|A|Ψ(µ) = ∂ λ E(λ, µ) λ=0 (2.14)
Thus, in this contribution, we use the same type of basis functions for A(µ) as for E(µ); see the results in Sec. 3.4. Note that computing Ψ(µ)|A|Ψ(µ) is not possible in DFT, without having a property-specific density functional [START_REF] Bauer | General operator ground-state expectation values in the Hohenberg-Kohn-Sham density-functional formalism[END_REF].

Numerical results

Guidelines

The quality of the corrections using Eqs. (2.7) and (2.11) is explored numerically. Technical details on the calculations are given in Appendix D.

The plots show the errors done by the approximations in the estimate of the energy: we choose a model, µ → E(µ), and let the empirical interpolation method choose which easier models (with weaker interactions) to extrapolate and get an approximation for E = E(µ = ∞). The plots show the error in the estimate of E made when considering approximations that use information only for μm ≤ µ. From the plots, we read off how small µ can be and still have "reasonable" accuracy. In thermochemistry, kcal mol -1 is a commonly considered unit, and is often considered as "chemical accuracy." For electronic excitations, one often uses eV units, and one often indicates it with one decimal place. "Chemical accuracy" is marked in the plots by horizontal dotted lines. The plots show the errors in the range of ±0.1 eV ≈ 2.3 kcal mol -1 .

We consider approximations using up to four points (thus chosen in [0, µ]). The first point µ 0 is always the value chosen µ 0 = µ shown on the xaxis of the plots, and the basis function associated to it is χ0 , the constant function; note that using only this pair (χ 0 , µ 0 ) corresponds to choosing E E(µ 0 )= the value provided by the model, i.e., no correction is applied. When the number of points is increased, further values of E(μ m ), chosen by the algorithm, are used with μm < µ.

The (maximal) parameter µ is considered between 0 and 3 bohr -1 . The model without correction (blue curve) reaches chemical accuracy for µ ≈ 3 bohr -1 for H -and harmonium, but only at µ ≈ 5 bohr -1 for H 2 in its ground state.

General behavior of errors

The plots in Fig. 2 for harmonium, H 2 , and H -have similar features and are discussed together. As the number of points used increases, the smallest value of µ for which the good accuracy is reached decreases. Note that FLEIM produces very small errors for values of µ larger than 2. However, with the chosen basis set, the algorithm presented in this contribution has difficulties correcting the errors for µ smaller than 1.

Possibility of error estimates

Some tests can be done to estimate the quality of the approximation. For example, we can compare how the approximations change when increasing the number of basis functions, K, in our approximation and consider |E K -E K-1 | as an asymptotically valid error estimate for E K-1 . One can notice in the above figures that, when the difference between, say, the 2- and the 3-point approximation error is larger than "chemical accuracy," so is the error in the 2-point approximation.

3.4.

Expectation values with FLEIM: r 12 and r 2 12 for harmonium

We look at the average distance between the electrons in harmonium. Figure 3(top) shows the error made by using Ψ(µ) instead of Ψ(µ = ∞) in computing the expectation value of r 12 , as well as the correction that can be achieved with FLEIM, using the same basis set as above (2.11). The inset in Figure 3(top) concentrates on the errors made in the region that could be considered chemically relevant (1 pm ≈ 0.02 bohr). We note the similarity with the behavior of in correcting E(µ).

Let us now examine the average square distance between the electrons, r 2 12 , in harmonium. While for computing the energy we explored correcting the missing short-range part of the interaction, we now ask whether it is possible to correct the error of using the model wave function, Ψ(µ) for the expectation value of an operator that is important at long range.

For ω = 1/2, we know the exact values of the expectation value of r 2 12 at µ = 0 and µ = ∞; they are 6 and (42

√ π + 64)/(5 √ π + 8) ≈ 8.21,
respectively (see, e.g., Ref. [START_REF] King | The electron correlation cusp[END_REF]). Note the large effect of the model wave function, Ψ(µ), in computing r 2 12 . Figure 3(bottom) shows the error made by using Ψ(µ) instead of Ψ(µ = ∞) in computing the expectation value of r 2 12 , as well as the effect of the correction that can be achieved with FLEIM, using the same basis set (2.11) as above. We note again the similarity with the behavior of in correcting E(µ) or the expectation value of the distance between electrons. The expectation value r 2 12 also illustrates another aspect: the effect of a change of the external potential on the energy. At first sight this may seem surprising, as the external potential is a one-particle operator, while r 2 12 is a two-particle operator. However, changing the one-particle operator also modifies the wave function and this affects the value of r 2 12 . In the case of harmonium, this can be shown analytically. Changing r 1 and r 2 to center-of-mass, R, and inter-particle distance, r 12 , cf. Appendix E, allows us to see explicitly that a modification of ω 2 , the parameter that specifies the external potential, affects the Schödinger equation in r 12 . It introduces a term proportional to ω 2 r 2 12 . The first order change in the energy when we change the external potential (ω 2 ) is thus proportional to r 2 12 . Our results in Fig. 3(bottom) show that our conclusions on model corrections are not modified by small changes in the external potential. Note that the centerof-mass Schrödinger equation also depends on ω 2 , but it is independent of µ and thus does not affect our discussion on model correction.

Comparison with DFAs

Instead of using extrapolation with FLEIM, one can use DFAs. While up to now the external potential did not change with µ, in DFA calculations a one particle potential that depends on µ is added in order to correct the density.

We consider here two DFAs, the local density approximation, LDA [START_REF] Savin | On degeneracy, near-degeneracy and density functional theory[END_REF][START_REF] Paziani | Local-spin-density functional for multideterminant density functional theory[END_REF], and one that reproduces that of Perdew, Burke and Ernzerhof (PBE) [START_REF] Goll | A short-range gradient-corrected density functional in long-range coupled-cluster calculations for rare gas dimers[END_REF][START_REF] Goll | Local and density fitting approximations within the shortrange/long-range hybrid scheme: application to large non-bonded complexes[END_REF] at µ = 0. Both approximations are modified to be µ-dependent. In particular, they vanish at µ = ∞.

As shown in Fig. 4(top) for harmonium, DFAs are clearly much better at small µ. However, they are not good enough. The figure suggests the range of µ for which DFAs are within chemical accuracy is similar to that obtained with the 3-point FLEIM. This is confirmed when comparing the results with DFAs and for the H 2 ; see Fig. 4(middle). Note that with FLEIM the errors at large µ are smaller.

Note also that the curves obtained with extrapolation are significantly flatter at large µ than those obtained with DFAs. This should not be surprising: DFAs transfer the large µ behavior, while extrapolation extracts it from information available for the system under study.

Furthermore, using ground-state DFAs for excited states does not only pose a problem of principle (questions its validity, as the Hohenberg-Kohn theorem is proven for the ground state), but can also show a deterioration of quality. However, there is no question of principle from the perspective of this contribution (of using a model and correcting it by extrapolation). Also, the error in the excited state seems comparable to that in the ground state, as seen in the example of the H 2 molecule, in the first excited state of the same symmetry as the ground state; see Fig. 4(bottom).

Conclusion and perspectives

In this contribution we have illustrated with a few models how to simplify the Hamiltonians by smoothly getting rid of the singularities in the system and thus have more numerically tractable problems. This simplification is obtained by introducing a parameter that, when it is equal to infinity, it corresponds to the original, plain Hamiltonian. After numerically solving a few simplified problems, the solution of interest is obtained by extrapolation.

We present a new (in the field) method for extrapolating the quantities of interest from few finite values (hence easy to solve) of the parameter by a technique borrowed from reduced basis paradigm: the empirical interpo-Étienne Polack, Yvon Maday, and Andreas Savin lation method. In contrast to DFAs, no parameters are fitted, no transfer from different system are made: only extrapolation is used. Note that in contrast to DFAs, improvement can be envisaged by either adding further points or using more appropriate basis functions and error estimates are asymptotically accessible.

Appendices A. On density functional approximations

In DFT, the existence of a universal functional of the density, F [ρ], i.e., the Hohenberg-Kohn theorem [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF], is rigorously proven [START_REF] Levy | Universal variational functionals of electron densities, firstorder density matrices, and natural spin-orbitals and solution of the vrepresentability problem[END_REF][START_REF] Lieb | Density functionals for Coulomb systems[END_REF]. However, obtaining accurately the value of F for a given density ρ(r ), while possible (see, e.g., Ref. [START_REF] Colonna | Correlation energies for some two-and fourelectron systems along the adiabatic connection in density functional theory[END_REF]), is exceedingly time-consuming. Computationally convenient DFAs exploit the knowledge of the density around a given point r in space. Typically,

F [ρ] ≈ R 3 f (ρ(r ), |∇ρ(r )|, . . . ) dr . (A.1)
The limitation of such an approach can be seen for a simple density functional, the Hartree term of the energy, *

E H [ρ] = R 6 ρ(r 1 )ρ(r 2 ) |r 1 -r 2 | (A.2)
when ρ(r ) = ρ A (r ) + ρ B (r ) and ρ A (r )ρ B (r ) ≈ 0 (i.e., ρ A and ρ B are spatially separated; their overlap decreases much faster than 1/|r 1 -r 2 |). Hohenberg and Kohn recognized the difficulty of approximating E H by expressions of the type given in Eq. (A.1) and suggested separating it from F [ρ]. However, this does not fundamentally solve the problem, as one can immediately see in one-electron systems, where E H has to be canceled by another term commonly expressed in DFAs by an ansatz of the forms given in Eq. (A.1). Note, that the problem would not exist for interactions that are not Coulomb (long-ranged), but short-ranged. For instance, if the interaction is Dirac's δ(r 1 -r 2 ) function, E H becomes exactly of the form of Eq. (A.1). For other short-range interactions one can approach such a form by using Taylor expansions. In recent years it has become popular to compensate for the limitation of the ansatz in Eq. (A.1) by adding "empirical" energy corrections to describe long-range effects. * The volume elements dr 1 dr 2 are omitted when the context is unambiguous.

Étienne Polack, Yvon Maday, and Andreas Savin

Another problem is that the antisymmetry of the electronic wave function is hidden in F [ρ]. As the most important effect of the antisymmetry is the Pauli repulsion, Kohn and Sham [START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF] proposed to consider the variational principle for a model in which particles do not interact. However, DFAs following the pattern of Eq. (A.1) are still in use. In fact, further separating terms from F [ρ] may even lead (for degenerate cases) to the question whether the limit of a noninteracting system is well-defined (see, e.g., Sec. 3.5 in [START_REF] Savin | On degeneracy, near-degeneracy and density functional theory[END_REF]).

B. Basis functions

In order to get an idea how the leading term of the correction behaves, we consider the missing part of the Hartree term,

Ψ(µ) W -W (µ) Ψ(µ) = 1 2 R 6 ρ(r 1 )ρ(r 2 ) 1 r 12 -w(r 12 , µ) . (B.1)
Separating the Hartree part from F [ρ] was already proposed by Hohenberg and Kohn [START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF], and it is also the dominant part in Ē(µ). Most molecular codes use Gaussian one-particle basis functions, so ρ is a linear combination of Gaussian functions. We consider a generic term, This integral is easily computed, e.g., by using Fourier transforms, and one obtains an expression that is proportional to the form of the basis function χj in Eq. (2.11), after arbitrarily relating 1/α 2 to the basis index j.

Let us now consider another basis set, constructed from the requirement that the functions decay as µ -2 and are finite at the origin. χ1,j (µ) = 1 + (j µ) 2 -1 (B.3)

where we choose again j = 1, 2, . . . , 10. The results are only slightly worse than those obtained when using the basis set of Eq. where we choose a j = 2 j/2 . The result for harmonium is shown in Fig. 5(bottom). We note a slight improvement of the results for the 2and 3-point approximations, as well as a change of sign of the error in the 4-point approximation. This direction of investigation deserves to be pursued.

C. The empirical interpolation method

The empirical interpolation method is a model-order reduction method introduced in [START_REF] Barrault | An 'empirical interpolation' method: Application to efficient reduced-basis discretization of partial differential equations[END_REF] as a way to efficiently find a reduced basis and approximate one particular function within a manifold of parameter dependent functions. The points at which to do the interpolation are referred as magic points [START_REF] Maday | A general multipurpose interpolation procedure: The magic points[END_REF]. For a family of basis functions χ i , i ∈ I with discrete points µ j , j ∈ J chosen on a regular grid close to zero (see Appendix D.2 for an analy-sis of the influence of the choice of the grid), we want to find a family of K functions and interpolation points with which to interpolate a test function f ∈ Span{χ i , i ∈ I}.

C.1. Algorithm for EIM

We assume that we have chosen some cost function C, e.g., a norm C[•] ≡ • , or, if we are, as in this contribution, only interested in correctly approximating the value for µ ≡ ∞, we choose the cost function as the absolute value at that extrapolation point

C[ϕ] = |ϕ(∞)| = lim µ→∞ |ϕ(µ)|.
First, select one of the basis functions. We can choose to add the constant function χ 0 , χ0 := χ 0 .

(C.1)

We then select the first interpolation point as the largest admissible µ available,

μ0 := max j∈J µ j . (C.2)
We can then define the first normalized interpolation function as

q 0 := χ0 χ0 (μ 0 ) . (C.3)
We can then create the first approximation with an interpolation scheme, for instance Lagrangian interpolation as

I 0 [f ] := f (μ 0 )q 0 . (C.4)
We now assume to have chosen K -1 functions χk , normalized functions q k . Let us also assume that we have selected K -1 interpolation points μk , k = 0, . . . , K -2. We define the I K-2 Lagrangian interpolation function as

I K-2 [f ] := K-2 k=0 β k q k , (C.5)
where the coefficients β k are determined by solving the system

K-2 k=0 β k q k (μ ) = f (μ ) for = 0, . . . , K -2. (C.6)
If we denote Ĩ as the set of indices of remaining basis functions and J as the set of indices of remaining interpolation points, we choose the next function as χK-1 := arg max

χ i , i∈ Ĩ C χ i -I K-2 [χ i ] , (C.7)
and the next interpolation point as μK-1 := arg max

µj , j∈ J χK-1 (µ j ) -I K-2 [ χK-1 ](µ j ) . (C.8)
We can then define the Kth normalised interpolation function as

q K-1 := χK-1 -I K-2 [ χK-1 ] χK-1 (μ K-1 ) -I K-2 [ χK-1 ](μ K ) . (C.9)
This system is represented by a lower triangular matrix with ones on the diagonal, and hence has a unique solution. The algorithm ends when the desired target accuracy is reached.

C.2. The forward looking empirical interpolation method (FLEIM)

To In this subsection, we compare on a simple function, the behaviour of EIM and FLEIM on a analytic test function E(µ) behaves like 1 + χ j (µ). The results for EIM are given in Fig. 6, and for FLEIM, in Fig. 7. First, we note that FLEIM provides better approximation together with more stable results when µ varies. Second, we note that the "wall" for j > 1 at small µ is also imporved. At large µ, all χ j have the same decay at large µ: this makes both methods fit to work in this regime.

D.2. Discretization for FLEIM

The interval between 0 and the value of µ under study was divided in 10 equal intervals (producing 11 points). FLEIM was used to select K(≤ 4) points and basis functions on which E(µ) was calculated.

We now investigate the effect of changing the grid of µ and use, for comparison a fixed finer grid of values of µ 1 , µ 2 , . . . ranging from 0 to the

D.4. Obtaining the model energy

In order to simplify the test of FLEIM, E(µ) was pre-calculated for a dense range of values µ and interpolated. The values for E(µ) were obtained with the program Molpro [START_REF] Amos | Molpro, a package of ab initio programs designed[END_REF] for H -and H 2 . This program was also used for the density functional calculations. For harmonium, it is possible to separate the variables in the Schrödinger equation. The center-of-mass equation can be solved exactly. The equation in r 12 was solved by discretization on a grid of 10 4 points between 0 and 10 bohr.

For H 2 the V5Z basis set of [START_REF] Dunning | Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen[END_REF] was used; for H -the aug-V5Z basis set of [START_REF] Woon | Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties[END_REF]. (The error in the energy of H -is too large if we do not augment the V5Z basis set with a diffuse basis function.) The aug-V5Z basis set was also used for the excited state of H 2 , as it has an important contribution of ionic states (H + • • • H -). For H 2 the equilibrium distance of 1.4 bohr was chosen for the ground state. For the excited state, the distance of 4.2 bohr was chosen. It is close to a minimum of the potential energy curve. Furthermore, this value can be compared with the accurate calculation of Ref. [START_REF] Nakashima | Solving the Schrödinger equation of hydrogen molecule with the free complement-local Schrödinger equation method: Potential energy curves of the ground and singly excited singlet and triplet states, Σ, Π, ∆, and Φ[END_REF].

Fig. 1 .

 1 Fig.1. Basis functions χj of Eq. (2.11), continuous curves with the color corresponding to j; and an (unknown) function to be approximated by linear combination on this basis (dot-dashed, gray). The unknown function in this figure is proportional to Ē(µ) of harmonium.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Errors for harmonium (top), H 2 (middle), and H -(bottom) using FLEIM with one to four points (1: blue curve, 2: brown curve, 3: green curve, 4: red curve). The abscissa represents the biggest µ allowed for use in the FLEIM algorithm. The error of the model without correction (blue curve) does not show up in the figure for the H 2 molecule because it is larger than the domain covered by the plot.

Fig. 4 .

 4 Fig.[START_REF] Flad | Explicit Green operators for quantum mechanical Hamiltonians. II. Edge-type singularities of the helium atom[END_REF]. Absolute errors for the harmonium molecule at equilibrium distance (top), for the H 2 molecule at equilibrium distance (middle), and in the first excited state of the same symmetry as the ground state (bottom): a µ-dependent LDA (black dashed curve), combined with the a µ-dependent Perdew-Burke-Ernzerhof approximation (PBE, gray dashed curve), FLEIM (3 points) (green curve), combined with a µ-dependent local density approximation. The abscissa represents the biggest µ allowed for use in the FLEIM algorithm. The insets zoom into the regions of small errors, the dotted line corresponding to the value of "chemical accuracy." Note the different ranges for ∆E.

( 2 .

 2 11), compare Fig. 2(top) and Fig. 5(top).Let us consider χ2,j (µ) = a j µ 1 + (a j µ) 3 -1 (B.4)

Fig. 5 .

 5 Fig.5. Errors for harmonium using FLEIM with up to four points (1: blue curve, 2: brown curve, 3: green curve, 4: red curve), using the basis sets given in (B.3) (top) and (B.4) (bottom). The abscissa represents the biggest µ allowed for use in the FLEIM algorithm.

  better adapt the method for extrapolation, we propose a double loop alternative: Instead of selecting sequentially first for a new basis function and then a new interpolation point -Eqs. (C.7) and (C.8) -we select the best pair ( χK-1 , μK-1 ) := arg max χ i , i∈ Ĩ arg min µj , j∈ J C χ i (µ j ) -I K-2 [χ i ](µ j ) . (C.10) D. Numerical details of the calculations D.1. Testing EIM and FLEIM with E(µ) = 1 + χ j (µ)

Fig. 6 .FLEIMFig. 7 .

 67 Fig. 6. Test of EIM using model E(µ) = 1 + χ j (µ).

Fig. 8 .

 8 Fig. 8. Results using a denser grid (lighter colors), for harmonium, H 2 at equilibrium distance, and H -; the curves obtained with the denser grid are shown in lighter colors. The error of the model without correction (blue curve) does not show up in the figure for the H 2 molecule because it is larger than the domain covered by the plot.

Funding Part of this work was supported by the French "Investissements d'Avenir" program, project ISITE-BFC (contract ANR15-IDEX-0003) ( ÉP). This work has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810367-project EMC2) (YM & AS).

Étienne Polack, Yvon Maday, and Andreas Savin

E. Change of coordinates in harmonium

Harmonium is characterized by the external potential given by Eq. (D.1). For two particles the variables r 1 , r 2 can be changed to those corresponding to the center of mass and the distance between particles,

yielding for the potential energy

The transformation of variables can be done also for the kinetic energy, and makes the model Schrödinger equation separable into a part (in R) that is independent of the model, and one (in r 12 ) that through 2.4 depends on the model (µ). Note that Ψ(µ)|r 2 12 |Ψ(µ) is also measuring the error due to the change of density with µ, for harmonium

Indeed, by using the Hellmann-Feynman theorem,

where R = (r 1 + r 2 )/2. As it is possible to separate variables R and r 12 in the Schrödinger equation,

Φ and φ both normalized to one, and

where ρ(r , µ) is the density of the model system, we have

Note that the change of r 2 12 with µ is only due to the change of the density with µ. Thus, the difference between 1 4 Ψ(µ)|r 2 12 |Ψ(µ) and 1 4 Ψ|r 2 12 |Ψ also indicates how much the density is affected by the model.


Van Long

Nguyen Huu

Julien Lallet
email: julien.lallet@nokia.com

Emmanuel Casseau

Laurent D'orazio
email: laurent.dorazio@irisa.fr

MASCARA-FPGA cooperation model: Query Trimming through accelerators

Keywords: FPGA acceleration, semantic caching, query optimisation

Processing sequences of online queries with an efficient response time is crucial for the new data management systems. To address this challenge, we proposed the ModulAr Semantic CAching fRAmework (MASCARA), in which Semantic Caching (SC) is used for a fast and accurate logs processing tool in a security monitoring system. Although MASCARA pointed out the expensive modules, we have not yet discussed an architecture in which they leverage the accelerators of Field Programmable Gate Array (FPGA). Therefore, in this paper, we propose a MASCARA-FPGA cooperation model and related major accelerators. We also explain the pipeline model in which multiple accelerators can run in parallel. Because Query Trimming can reduce the response time by factors of up to 2.8x for a single kernel, we confirm the potential of building the full MASCARA-FPGA model.

INTRODUCTION

In recent years, Big Data has led to the emergence of new data management systems as an alternative to Relational Database Management Systems (RDBMS), such as: column-stores, or Massive Parallel Processing (MPP). Some of these systems can be deployed on a specified query engine, such as Spark [START_REF] Armbrust | Spark SQL: Relational Data Processing in Spark[END_REF], on top of a distributed file system. As a case study, an IoT security monitoring system with the HTTP log analysis tool is designed to ensure the safety by observing unusual events. This monitoring tool is based on a variety of specified information reasoning from an enormous number of access logs. For example, this tool can consist of select-project queries to get the traffic for a given time. Thus, reducing the response time of these online queries is a real challenge.

To achieve this goal, there are two basic approaches: 1) Speeding up the execution time of a query by offloading it from CPU to a high-speed computing unit; or 2) Exploiting resources and knowledge contained in the queries themselves by a caching system. With the first approach, to overcome the CPU limits, the trend is evolving to explore the acceleration of query execution via Field Programmable Gate Arrays (FPGAs) [START_REF] Fang | In-memory database acceleration on FPGAs: a survey[END_REF]. FPGAs with partial runtime-reconfiguration capability can close the gap between flexibility of CPU and the performance of specialized hardware accelerators. In the mean time, with the second approach, the query engines have not yet implemented any technique to exploit the semantic relevance between 𝑄1 and 𝑄2. Therefore, our goal is to build a query processing platform on FPGA with semantic caching (SC) as a combination of these two approaches to optimize the performance.

Related works On one hand, we have some noticeable accelerating approaches on FPGA without constructing a cache or applying the semantic theory. Sukhwani et al. [START_REF] Sukhwani | A Hardware/Software Approach for Database Query Acceleration with FPGAs[END_REF] proposed FGPA accelerated database system in which the analytics queries are offloaded from the host running the workload to the FPGA composable accelerators. Mueller et al. [START_REF] Mueller | Glacier: A Query-to-Hardware Compiler[END_REF] proposed Glacier, the query to hardware compiler, which can compose a digital circuit using various implemented operators. This approach is only suitable for a known query set. Dennl et al. [START_REF] Dennl | Acceleration of SQL Restrictions and Aggregations through FPGA-Based Dynamic Partial Reconfiguration[END_REF] presented concepts for multi-processing queries in MySQL. This approach focuses on restriction and aggregation operators and thus cannot evaluate more complex queries. Sidler et al. [START_REF] Sidler | DoppioDB: A Hardware Accelerated Database[END_REF] presented DoppioDB, which extends MonetDB on FPGAs. In cooperation with DoppioDB, Owaida et al. [START_REF] Owaida | Centaur: A Framework for Hybrid CPU-FPGA Databases[END_REF] provided the Centaur framework to bridge the gap between CPUs and FPGAs. Recently, Wang et al. [START_REF] Wang | Relational query processing on OpenCL-based FPGAs[END_REF] approached the advantage of programmability offered by the OpenCL to generate an optimal plan for relational query processing on FPGA.

On the other hand, there are some remarkable researches or proposals regarding semantic exploitation. SC was studied and proved by Dar et al. [START_REF] Dar | Semantic Data Caching and Replacement[END_REF] or Keller et al. [START_REF] Keller | A Predicate-Based Caching Scheme for Client-Server Database Architectures[END_REF] as an effective solution compared to pages or tuples caching. Moreover, we have the research of Zhou et al. [START_REF] Zhou | Automated Verification of Query Equivalence Using Satisfiability modulo Theories[END_REF] to exploit the semantic equivalence between queries. Although these researches were not presented in the context of acceleration on FPGA, they confirmed that data reasoning based on the semantic approach could speed up the performance of the data management system. Meanwhile, Salami et al. [START_REF] Salami | Accelerating Hash-Based Query Processing Operations on FPGAs by a Hash Table Caching Technique[END_REF] proposed an efficient cache on FPGA to store partly or entirely the hash table for hash join operator where CPU is difficult to have strong robustness and high speed. Unfortunately, their cache did not exploit the advantages of SC or query equivalence.

Our previous works. To the best of our knowledge, proposing SC on FPGA has not yet been studied. Thus, we had a vision of multi-layer framework for IoT security monitoring tool [START_REF] Orazio | Semantic Caching Framework: An FPGA-Based Application for IoT Security Monitoring[END_REF] without addressing the semantic processing. Later, we presented ModulAr Semantic CAching fRAmework (MASCARA) [START_REF] Van Long Nguyen Huu | MASCARA (ModulAr Semantic CAching fRAmework) towards FPGA Acceleration for IoT Security Monitoring[END_REF] that describe the semantic processing on host only. We also pointed out which module in MASCARA is expensive in computing time so that it can be interesting to offload and accelerate by FPGA.

Contributions. This paper address the following challenge: we have not yet had an MASCARA architecture from the hardware perspective in which the modules are expressed on FPGA to leverage its computing kernels. Therefore, in order to address this challenge, our objective is to present the following new contributions:

(i) A cooperation model in which MASCARA's stages or modules are addressed on FPGA as computing kernels.

(ii) Query Trimming kernels in the model of pipeline execution .

(iii) The preliminary results from single to multiple kernels to confirm the acceleration of Query Trimming on FPGA.

MASCARA ARCHITECTURE

This section presents briefly the SC as an optimization solution in which the modules of MASCARA are described.

Semantic Caching. SC enables data reasoning effectively, reducing data transmission on servers, and exploiting the relevance of query [START_REF] Dar | Semantic Data Caching and Replacement[END_REF], [START_REF] Keller | A Predicate-Based Caching Scheme for Client-Server Database Architectures[END_REF]. Consider that SC stores a semantic entity SE: (𝑑𝑎𝑡𝑒 < ′ 17/03/2020-10 : 00 : 00 ′ ) that points to a result set of data in cache. Then, an incoming query INCQ: (𝑖𝑝 = ′ 208.115.111.72 ′ ) would be divided into: PROBEQ: (𝑑𝑎𝑡𝑒 < ′ 17/03/2020 -10 : 00 : 00 ′ ∧𝑖𝑝 = ′ 208.115.111.72 ′ ) processed in cache to retrieve the portion of the result available. Meanwhile, REMAINQ: (𝑑𝑎𝑡𝑒 ≠ ′ 17/03/2020 -10 : 00 : 00 ′ ∧ 𝑖𝑝 = ′ 208.115.111.72 ′ ) has to retrieve any missing tuples from the server.

MASCARA overview. Based on this approach, the MASCARA [START_REF] Van Long Nguyen Huu | MASCARA (ModulAr Semantic CAching fRAmework) towards FPGA Acceleration for IoT Security Monitoring[END_REF] is composed of 4 main stages: (1) Query Broking, (2) Query Trimming, (3) Semantic Management, and (4) Result refining. The modular approach allows us to reuse cache services or focus only on developing a cache solution. In MASCARA's previous work, all of the modules are deployed on the host at the software level.

Within MASCARA, the Query Broking stage materializes the query into various types of elements. The elements of an incoming query would be turned and stored as a tuple 𝑄 < 𝑄 𝑅 , 𝑄 𝐴 , 𝑄 𝑃 , 𝑄 𝐹 , 𝑄 𝐶 > with R is a Relation, A is an Attribute, P is a Predicate, F is an aggregate Function (SUM, MAX, MIN), and C is the number of records. This semantic tuple 𝑄 will be trimmed to Probe Queries PQ(s) and Remainder Queries RQ(s) by Query Trimming stage based on its relationship with semantic segments 𝑆 < 𝑆 𝑅 , 𝑆 𝐴 , 𝑆 𝑃 , 𝑆 𝐹 , 𝑆 𝐶 > that are stored in Semantic Management stage. According to SC theory, each 𝑆 𝑖 points to the result set of data that is presented by Semantic Data Region in the cache. Last, the Result Refining stage will create the final answer for the user by combining the result sets of PQ(s) and (RQ)s together. Query Trimming. This stage plays an important role due to its complex procedure of rewriting the incoming query to PQ(s) and RQ(s) [START_REF] Van Long Nguyen Huu | MASCARA (ModulAr Semantic CAching fRAmework) towards FPGA Acceleration for IoT Security Monitoring[END_REF]. Thus, we proposed two sub-stages: Semantic Matching and Semantic Extracting, that could be executed in parallel.

Semantic Matching is a process in which we have Attribute Matching of (𝑄 𝐴 , 𝑆 𝐴 ) and Predicate Matching of (𝑄 𝑃 , 𝑆 𝑃 ) from the input pair (𝑄, 𝑆). The Attribute Matching computes the Common Attribute (𝐶 𝐴 = 𝑄 𝐴 ∩ 𝑆 𝐴 ) and Difference Attribute (𝐷 𝐴 = 𝑄 𝐴 -𝑆 𝐴 ) from (𝑄 𝐴 , 𝑆 𝐴 ). Meanwhile, the Predicate Matching identifies the Intersection (𝑄 𝑃 ∩ 𝑆 𝑃 ), Difference (𝑄 𝑃 \ 𝑆 𝑃 ) and Implication (𝑄 𝑃 → 𝑆 𝑃 ) relationship from (𝑄 𝑃 , 𝑆 𝑃 ).

Based on the results of them, the output of Semantic Matching could be one of the following case: Totally Contained, Horizontal Matching, Vertical Matching, Mixed Matching and NO Matching that are presented as semantic box approach in [START_REF] Dar | Semantic Data Caching and Replacement[END_REF]. For each case, the Semantic Extracting will create the PQ(s) and RQ(s) as outputs.

Limitation. Up to now, MASCARA hasn't addressed the stages or modules as accelerated kernels on FPGA, such as: the Query Trimming engine, or the cache in Semantic Management. Thus, to address this challenge, our vision is to propose a cooperation model in which we leverage the computing kernels of FPGA.

MASCARA-FPGA COOPERATION MODEL

In this section, we introduce a hardware perspective model to address the cooperation between MASCARA and relevant major components on FPGA for semantic processing and query executing (as shown in Figure 2). Then, we present the accelerators of Query Trimming within a pipeline execution model. Query Reformer. Besides the MASCARA's stages, we also have Query Reformer engine to remove the unnecessary elements in queries before materializing them into 𝑄 by Query Broking stage. As a result, by improving the appearance of the query, this engine will reduce the complexity of generated remainder queries later.

FPGA-Adapter. To maintain the communication between modules on host and FPGA, we use an adapter to transfer 𝑄 and also the data blocks through the PCIe interface. This FPGA-Adapter has to compose a foreign-data wrapper, i.e., Java Wrapper, to build the calling method of accelerators. MASCARA on FPGA. Because the Query Trimming is expensive in computing time compared to the other stages of MASCARA, we want to offload it to FPGA [START_REF] Van Long Nguyen Huu | MASCARA (ModulAr Semantic CAching fRAmework) towards FPGA Acceleration for IoT Security Monitoring[END_REF]. Moreover, Query Trimming needs to cooperate with Semantic Management to generate the output PQ(s) and RQ(S). Thus, placing also the Semantic Management on FPGA is a complementary solution to reduce the data communication between host and FPGA. Our vision is implementing the modules of Query Trimming and Semantic Management on FPGA using the same concepts of query rewriting, cache organization and replacement strategy that are presented in our previous work [START_REF] Van Long Nguyen Huu | MASCARA (ModulAr Semantic CAching fRAmework) towards FPGA Acceleration for IoT Security Monitoring[END_REF].

For Query Trimming, we construct different sets of accelerators that are relevant to semantic matching and extracting, such as accelerators of Attribute Matching or Predicate Matching. They can be invoked from the host or other accelerators directly according to execution data flow. Meanwhile, we place the Semantic Management on DDR of FPGA and divide it into two parts: Semantic Info Table to store list of 𝑆 𝑖 and Semantic Data Region to store relevant result set of 𝑆 𝑖 . These DDRs have a high bandwidth to and from the configurable region of MASCARA to ensure the performance of Query Trimming's kernels. For each 𝑆 𝑖 , we support to store the element < 𝑆 𝑃 > in forms of Disjunctive Normal Form (DNF) in which we have 𝑃 1 ∨ 𝑃 2 ∨ ... ∨ 𝑃 𝑛 . In detail, a predicate 𝑃 𝑖 is separated to the Range Indicators 𝑅𝐼 in forms of Conjunctive Normal Form (CNF): 𝑃 𝑖 = 𝑅𝐼 1 ∧ 𝑅𝐼 2 ∧ ... ∧ 𝑅𝐼 𝑛 . A 𝑅𝐼 is an object that could be represented as an interval (𝐿𝑜𝑤𝑒𝑟 𝐵𝑜𝑢𝑛𝑑, 𝑈 𝑝𝑝𝑒𝑟𝐵𝑜𝑢𝑛𝑑) or an 𝐸𝑥𝑎𝑐𝑡𝑉 𝑎𝑙𝑢𝑒. Currently, we provide the comparison operators: <, >, =, ≥, ≤, ≠ for RI.

Query Process Controller. QPC controls the trimming life cycle of 𝑄 from Semantic Matching to Semantic Extracting accelerators. To leverage the most out of FPGA's power, PQ(s) and RQ(s) are considered to be accelerated by using the library of database (DB) operators within the QPC's management. If PQ(s) or RQ(s) contains any unsupported operator, they would be transferred and executed on the host instead. Of course, for the execution of RQ(s), we take into account the challenge of traditional data offloading mechanism where the overhead of data transmission from host to FPGA accelerators can occur.

Based on the status of Operator Manager, the Finite State Machine (FSM) creates the control signals to the Interconnection Unit in which a physical planning is constructed for query executing.

MASCARA accelerators

We present the data flow of Query Trimming on FPGA in Figure 3 to emphasize the effectiveness of accelerators, such as: Attribute Matching, Predicate Matching, and Semantic Extracting. Thus, the most advantage of FPGA kernels is that the throughput of Query Trimming (PQ, RQ) increases by processing each pair (𝑄, 𝑆 𝑖 ) in parallel from matching to extracting.

Attribute Matching kernel. We take 𝑄 < 𝑄 𝐴1,𝐴2,...,𝐴𝑛 > from 𝑄 and 𝑆 < 𝑆 𝐴1,𝐴2,...,𝐴𝑚 > from 𝑆 𝑖 as two inputs. Using a nested iteration, we compare the attributes of two input arrays. As a result, the outputs could be maintained as the collection: 𝐶 𝐴 and 𝐷 𝐴 .

Predicate Matching kernel. We take 𝑄 < 𝑄 𝑃 1,𝑃 2,...,𝑃𝑛 > and 𝑆 < 𝑆 𝑃 1,𝑃 2,...,𝑃𝑚 > as two inputs. The Predicate Matching accelerator is more complex than the Attribute Matching accelerator. We have to implement three engines: Intersection, Difference, and Implication from predicates to ranges. Predicate Matching is expensive due to the complexity of Difference engine when its computing becomes intense quickly with a large number of 𝑄 𝑃 or 𝑆 𝑃 .

Semantic Extracting kernel. Using the outputs from Attribute and Predicate Matching, we can construct PQ and RQ(s) in the Semantic Extracting kernel. For example, 𝐷 𝐴 is used to build 𝑅𝑄1 < 𝑅𝑄 𝐴 > meanwhile 𝐷 𝑃 is used to build 𝑅𝑄1 < 𝑅𝑄 𝑃 >. As soon as PQ and RQ(s) appear, they will be stored in the buffer of execution and wait for the decision of QPC that they could be executed/accelerated on FPGA or not based on the available DB accelerators.

Pipeline model for multiple accelerators

Because these MASCARA accelerators can be executed in parallel, we gain two main benefits: 1) We can apply a pipeline model to process the pair of (𝑄, 𝑆 𝑖 ) as soon as the accelerators are available; 2) Suppose that hardware resources of FPGA are sufficient to create the multiple instances of accelerators, we could reduce the total matching time between 𝑄 and a set of 𝑆 𝑖 .

Assuming that we have two instances for each MASCARA's accelerators. From a theoretical point of view, we consider the execution times of accelerators are equal without regarding the complexity of their prototypes. As shown in Figure 4, we can execute (𝑄1, 𝑆 1 ), (𝑄1, 𝑆 2 ) in parallel and forward result to Extracting as next stage of pipeline. Thus, the throughput of Query Trimming with multiple accelerators in the pipeline model will be increased compared to the sequential or purely pipeline version. Meanwhile, from the practical point of view, to build this pipeline model for multiple accelerators, we have to cope with defining the streaming connections and managing the data flow between pipeline stages.

PRELIMINARY EXPERIMENTS

The completely MASCARA-FPGA becomes a real challenge in which we need the kernels to accelerate: 1) Query Trimming, and 2) Semantic Data Management for query execution. Therefore, to show the potential of our model, in this experiment, we currently focus on confirming the capability of Query Trimming on FPGA without accelerating their outputs (PQs, RQs) based on Semantic Data Management.

Configuration. The host is an Intel® Core™ i7-6700 CPU 3.40GHz, and FPGA is an Alveo U200 Xilinx Card. To develop the accelerators, we use the Vitis™ platform from Xilinx.

Workloads. Because we have not accelerated yet the PQ(s) or RQ(s) from Query Trimming, we only need to construct a list of segments 𝑆 to match and extract with query 𝑄. Thus, we initialize Semantic Info Table on DDR with 100 segments that has up to 15 attributes based on EDGAR log file data set of SEC.gov [START_REF]Securities and exchange commission[END_REF]. We create Evaluation. An advantage of Query Trimming is that we do not need too much hardware resources. In detail, it consumes 5703 FFs, 7535 LUTs and 31 BRAMs. We compare the response time of Query Trimming through host only version and multiple kernels on FPGA. As shown in Figure 5, the host only version (green) always has the lowest performance. In particular, if we just use only one Query Trimming kernel (blue), the average response time is reduced by factors of up to 2.8x. Moreover, we confirm that the more kernels we have, the better speed-up we can get. For example, with 4 kernels (red), we have a better response time by factors of up to 9.6x.

The management of multiple kernels was presented through Figure 4 in section 3.2. In complete MASCARA-FPGA later, although we can create more kernels, we have to optimize their instances in context of the outputs execution (PQs, RQs).

Fig. 1 .

 1 Fig. 1. An overview of MASCARA architecture

Fig. 2 .

 2 Fig. 2. MASCARA-FPGA co-operation model and its major components. The stages of MASCARA are expressed in red font, i.e. Semantic Matching Accelerator

Fig. 3 .

 3 Fig. 3. Dataflow of Query Trimming on FPGA

Fig. 4 .

 4 Fig. 4. Query Trimming kernels in pipeline model

Fig. 5 .

 5 Fig. 5. Total response time of Query Trimming stage

CONCLUSION & FUTURE WORKSIn this paper, we propose the MASCARA-FPGA cooperation model in which stages or modules are addressed as computing kernels. Through preliminary results of Query Trimming, we confirm that building a completely MASCARA-FPGA is potential.Future work includes the full implementation of MASCARA-FPGA for Select-Project-Join queries with Aggregation. Moreover, we have to study the effectiveness of data organization, i.e, dictionary-based approach instead of raw data, or replacement strategy for Semantic Data Region on FPGA.