Sabrina Gaba
email: sabrina.gaba@dijon.inra.fr

François Bretagnolle

Thierry Rigaud

Laurent Philippot

Vesna Gagic

Aidan M Keith

Managing biotic interactions for ecological intensification of agroecosystems

Keywords: agroecology, pest control, land use, adaptation, ecosystem services, biodiversity

Agriculture faces the challenge of increasing food production while simultaneously reducing the use of inputs and delivering other ecosystem services. Ecological intensification of agriculture is a paradigm shift, which has recently been proposed to meet such challenges through the manipulation of biotic interactions. While this approach opens up new possibilities, there are many constraints related to the complexity of agroecosystems that make it difficult to implement. Future advances, which are essential to guide agricultural policy, require an eco-evolutionary framework to ensure that ecological intensification is beneficial in the long term.

INTRODUCTION

The manipulation and regulation of biotic interactions and the functioning of agroecosystems in order to increase agricultural production sustainably while maintaining the diverse services provided by these ecosystems is the major challenge of the Twenty first century. Agricultural lands are more complex and diverse ecosystems than pictured in the classical cliché. Most of this biodiversity is hidden in the soil, with up to 10 7 prokaryotes per gram of soil and at least 10 5 bacterial species [START_REF] Torsvik | High diversity in DNA of soil bacteria[END_REF][START_REF] Whitman | Prokaryotes: the unseen majority[END_REF][START_REF] Gans | Checklist and bibliography of the Amphisbaenia of the world[END_REF], and agroecosystems also host insects, weeds, birds and rodents, leading to complex interactions and trophic networks (Figure 1). Managing interactions between and within the different trophic levels may enhance agroecosystem functions and increase crop productivity and improve pest control and sustainability [START_REF] Mcneely | Ecoagriculture, Strategies to Feed the World and Save Wild Biodiversity[END_REF][START_REF] Bommarco | Ecological intensification: harnessing ecosystem services for food security[END_REF]Tixier et al., 2013a). For example, interactions between belowground soil organisms and aboveground plant communities influence ecosystem processes and properties with consequences for the provision of the related services [START_REF] Bardgett | A temporal approach to linking aboveground and belowground ecology[END_REF]. The direct management of biotic interactions such as the use of microbial strains to promote plant growth [START_REF] Shennan | Biotic interactions, ecological knowledge and agriculture[END_REF][START_REF] Doré | Facing up to the paradigm of ecological intensification in agronomy: revisiting methods, concepts and knowledge[END_REF][START_REF] Ekstrom | Pest control in agro-ecosystems: an ecological approach[END_REF] and the enhancement of services that ecosystems provide naturally, such as pest control, underlies the ecological intensification of cropping systems (Figure 2). Given its potential as an alternative to chemical intensification, interest in ecological intensification has been growing, both in ecological and agronomic research [START_REF] Cassman | Ecological intensification of cereal production systems: yield potential, soil quality, and precision agriculture[END_REF][START_REF] Bommarco | Ecological intensification: harnessing ecosystem services for food security[END_REF][START_REF] Garnett | Sustainable intensification in agriculture: premises and policies[END_REF]. However, current research often relies on the manipulation of interactions between a small number of species, for example a pest and its predator, even though these interactions are part of a complex network. This overlooks (i) the possible resulting trophic cascades within the network, (ii) the possible eco-evolutionary responses of the species in the network, and (iii) the resistance and resilience of the ecosystem.

In this review, we argue that, whatever the methods used for ecological intensification (e.g., direct management of biotic interactions or intensification of the biological processes supporting ecosystem services), there is a need for a deeper understanding of biological interactions within agroecosystems to avoid techniques that are ineffective or, worse, have negative impacts that outweigh the benefits to the services. Such knowledge would gain to be investigated within an eco-evolutionary framework, to explore effects in the medium term (effectiveness of the management techniques) and the long term (sustainability of the management techniques). This would also allow taking advantage of the possibilities provided by complex ecological modeling (e.g., network models in the widest sense from food webs to landscape networks which take account of complex ecological properties to test management strategies, Tixier et al., 2013b). The management of the multiple services delivered by agroecosystems and of the corresponding trade-offs has recently been addressed by [START_REF] Bommarco | Ecological intensification: harnessing ecosystem services for food security[END_REF], while [START_REF] Tscharntke | Global food security, biodiversity conservation and the future of agricultural intensification[END_REF] outlined the shortcomings of a dichotomic view between land sharing and land sparing for real world application. In this review, we focus on how to ensure the effectiveness and sustainability of ecological intensification. After introducing the different definitions of ecological intensification, we present significant achievements in the field. We then suggest a conceptual framework for successful management of biotic interactions. Finally, we discuss approaches and tools for practical applications of ecological intensification.

DEFINITIONS OF ECOLOGICAL INTENSIFICATION

The term "ecological intensification" is relatively recent and there is still no generally accepted meaning [START_REF] Garnett | Sustainable intensification in agriculture: premises and policies[END_REF]. This shows that agricultural farmlands provide a suitable habitat for species other than the main crop, opening up new possibilities for resource and pest control by using the natural functions of farmland biodiversity.

"Ecological intensification" was originally defined as the process of increasing crop production to satisfy future food demand while meeting acceptable standards of environmental quality [START_REF] Cassman | Ecological intensification of cereal production systems: yield potential, soil quality, and precision agriculture[END_REF][START_REF] Cassman | Scientific challenges underpinning the food-versus-fuel debate[END_REF][START_REF] Garnett | Sustainable intensification in agriculture: premises and policies[END_REF]. Its aim is to maximize the primary production per unit area without compromising the ability of the system to sustain its productive capacity. However, ecological intensification can be more than simply a form of agricultural intensification that preserves the environment by making best use of environmental goods and services [START_REF] Pretty | Agricultural sustainability: concepts, principles and evidence[END_REF][START_REF] Bonny | Ecologically intensive agriculture: nature and challenges[END_REF]. It is based on management of ecosystem processes rather than fossil fuel inputs by (i) integrating biological and ecological processes into food production processes and (ii) minimizing the use of non-renewable inputs. In this sense, "ecological intensification" of agriculture is based on the principles of persistence and resilience of the agroecosystem and the principle of autarchy, i.e., the capacity to deliver outputs from inputs and resources acquired from within the system boundaries [START_REF] Pretty | Agricultural sustainability: concepts, principles and evidence[END_REF][START_REF] Power | Ecosystem services and agriculture: tradeoffs and synergies[END_REF]. One way to maximize production through ecological intensification is to manage organisms providing regulation services such as pest control, pollination and soil nutrients [START_REF] Bommarco | Ecological intensification: harnessing ecosystem services for food security[END_REF].

More recently, "ecological intensification" has been placed in the broader concept of ecological engineering which is the use of ecological science and theory to design, construct and manage ecosystems to benefit both humanity and nature [START_REF] Mitsch | Ecological engineering: a field whose time has come[END_REF]. One of the key principles of this paradigm is to integrate biological and ecological processes into food production, focusing on the agricultural system and landscape design. It relies on the premise that natural ecosystems may provide good models for designing agricultural systems that match environmental, social and sustainability objectives.

MANAGING BIOTIC INTERACTIONS: SUCCESSES AND POTENTIAL

Ecological intensification aims to promote beneficial biological interactions to limit the massive use of chemical inputs such as pesticides and fertilizers and to reduce the environmental impact. A conventional approach for decreasing pesticide use is based on antagonistic interactions between pests and their natural enemies. One of the flagship successes of such biological control in crops was the release of the parasitoïd Encarsia formosa against the greenhouse whitefly Trialeurodes vaporariorum in the 1970s as a substitute for pesticides to which the whitefly had evolved to be resistant [START_REF] Hoddle | Biology and use of the whitefly parasitoid Encarsia formosa[END_REF][START_REF] Bale | Biological control and sustainable food production[END_REF]. The importance of antagonistic interactions in crop fields is illustrated by the so called "disease-suppressive soils" in which crop plants suffer less from certain diseases even though the virulent soil-borne pathogens are present. In these soils, interactions between the pathogens and the indigenous soil microorganisms are responsible for reducing the severity of the disease by inhibiting the growth or activity of the pathogen [START_REF] Mendes | Deciphering the rhizosphere Microbiome for disease-suppressive bacteria[END_REF]. The potential for biocontrol of diseases by using interactions such as amensalism (e.g., antibiosis), competition and parasitism between fungal or bacterial strains and plant pathogens has been widely studied in microbial ecology [START_REF] Haas | Biological control of soil-borne pathogens by fluorescent pseudomonads[END_REF][START_REF] Harman | Overview of mechanisms and uses of Trichoderma spp[END_REF]Philippot et al., 2013a).

Similarly, positive interactions such as facilitation and mutualism can also help to reduce fertilizer inputs and increase primary production. In an intercropped agroecosystem, the exudation of phosphorous-mobilizing compounds or other mechanisms for increasing phosphate availability by one species can facilitate phosphorous uptake by the other species [START_REF] Hinsinger | P for two, sharing a scarce resource: soil phosphorus acquisition in the rhizosphere of intercropped species[END_REF]. Saprophytic microorganisms thriving in the rhizosphere also increase soil resource availability by cycling nutrients, which is beneficial for the plants. Plant growth promoting rhizobacteria (PGPR), which are commonly applied in developing countries where the use of mineral fertilizer is limited by costs, can also trigger growth by synthesizing volatiles mimicking plant hormones [START_REF] Ryu | Bacterial volatiles promote growth in Arabidopsis[END_REF]. Managing plant-microbe interactions is of interest not only for optimizing nutrient cycling but also for reducing consequences such as leaching and greenhouse gas emissions. It has recently been demonstrated that compounds inhibiting the microbial oxidation of ammonium into nitrate released by plant roots resulted in reduced field emissions of N 2 O [START_REF] Subbarao | Evidence for biological nitrification inhibition in Brachiaria pastures[END_REF]. Selection for such traits may lead to a new generation of crop cultivars that limit microbial emissions of greenhouse gases [START_REF] Philippot | Towards food, feed and energy crops mitigating climate change[END_REF].

PROMOTING BIODIVERSITY TO ENSURE SUSTAINABILITY

Promoting biodiversity in agroecosystems can increase their sustainability in various ways. Firstly, it is now evident that biodiversity is important for ecosystem functioning and this has been demonstrated for several communities [START_REF] Isbell | High plant diversity is needed to maintain ecosystem services[END_REF][START_REF] Cardinale | Biodiversity loss and its impact on humanity[END_REF]. In an early study, [START_REF] Tilman | Productivity and sustainability influenced by biodiversity in grassland ecosystems[END_REF] showed that plant productivity increased with plant species richness. Similarly, predator diversity strengthens pest control of either single or multiple species of prey [START_REF] Snyder | Predator biodiversity strengthens aphid suppression across single-and multiple-species prey communities[END_REF][START_REF] Takizawa | Predator biodiversity increases the survivorship of juvenile predators[END_REF] and plant diversity can have beneficial effects on pest control by encouraging the natural enemies of crop pests, either because non-crop plants increase the number and diversity of enemies [START_REF] Perdikis | Enhancing pest control in annual crops by conservation of predatory Heteroptera[END_REF] or because they increase the number of links in the food (or interaction) webs. In this second case, an increased number of links increases the stability of networks including natural enemies, thus reducing the risk of enemy extinction [START_REF] Macfadyen | Landscape structure influences modularity patterns in farm food webs: consequences for pest control[END_REF]. Diversifying cropping systems with non-crop species may improve top-down control of natural herbivores. In a recent paper, [START_REF] Balmer | Noncrop flowering plants restore top-down herbivore control in agricultural fields[END_REF] found that non-crop flowering plants planted within or around Brassica oleracea fields had a significant effect on the complex plant-herbivore-parasitoid-predator food web through the increased rate of parasitism and predation on the target cabbage pest species. Plant production can also be affected by microbial community composition and diversity [START_REF] Wardle | Ecological linkages between aboveground and belowground biota[END_REF][START_REF] Maherali | Influence of Phylogeny on fungal community assembly and ecosystem functioning[END_REF]. It was also demonstrated experimentally that microbial diversity loss was detrimental to the cycling of nitrogen, a nutrient essential for plant growth [START_REF] Philippot | Loss in microbial diversity affects nitrogen cycling in soil[END_REF], and that rare soil microbes may play a role in crop protection by improving aboveground and belowground plant defenses [START_REF] Hol | Reduction of rare soil microbes modifies plant-herbivore interactions[END_REF]. Secondly, it is argued that increased biodiversity can also "insure" ecosystems against decline in functioning under fluctuating conditions, thus increasing ecosystem stability [START_REF] Yachi | Biodiversity and ecosystem productivity in a fluctuating environment: the insurance hypothesis[END_REF]. This "insurance hypothesis" is based on the intuitive concept that the probability of finding species able to adapt to changing conditions and allowing ecosystem functioning is greater in a more diverse ecosystem.

Decreased soil biodiversity reduces the resilience of organic matter decomposition to a disturbance caused by adding copper [START_REF] Griffiths | Ecosystem response of pasture soil communities to fumigationinduced microbial diversity reductions: an examination of the biodiversityecosystem function relationship[END_REF]. Furthermore, manipulation of both richness and evenness of a microbial guild involved in nitrogen cycling showed that its resistance to salinity disturbance was lower when the initial communities were very uneven or dominated by a few species, suggesting that evenness can also be important for ecosystem stability [START_REF] Wittebolle | Initial community evenness favours functionality under selective stress[END_REF]. Finally, it was suggested that the influence of biodiversity on ecosystem functioning and stability depends on the strength of the interactions between species [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Rooney | Integrating diversity, food web structure and stability[END_REF]. Compartmentalized trophic networks, which consist of groups of species that have a higher probability of interacting with each other rather than with species of other groups, showed increased stability against perturbation [START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF]. The type of interaction is also important since, in contrast to trophic networks, stability is increased by a densely connected, nested structure in mutualistic networks [START_REF] Thebault | Stability of ecological communities and the architecture of mutualistic and trophic networks[END_REF].

While increasing biodiversity is a promising way of enhancing biotic interactions to reduce the use of chemical inputs, the trade-offs and synergies of these biotic interactions are rarely quantified. For example, complementarity between predators and niche partitioning is often suggested as an explanation of the higher success of a diversity of predators on resource consumption (e.g., [START_REF] Northfield | Niche saturation reveals resource partitioning among consumers[END_REF] and, by extension, pest control. However, increasing the diversity of predators may increase the frequency of intraguild predator species and, therefore, diminish the effectiveness of pest control through predator-predator interactions (e.g., [START_REF] Schausberger | Combined versus single species release of predaceous mites: predator-predator interactions and pest suppression[END_REF]. Furthermore, although increasing plant diversity can be translated into better pest control, it is often less effective than intensive chemical treatment (e.g., [START_REF] Brown | Role of biodiversity in integrated fruit production in eastern North American orchards[END_REF], leading to a significant decrease in the overall crop yields [START_REF] Letourneau | Does plant diversity benefit agroecosystems? A synthetic review[END_REF]. Therefore, a better knowledge of the network properties is still needed to understand the relationship between complexity and ecological stability [START_REF] Montoya | Ecological networks and their fragility[END_REF] and to make efficient use of increased biodiversity as an ecological service.

IMPORTANCE OF THE SPATIAL SCALE AND LANDSCAPE ORGANIZATION FOR BIOTIC INTERACTION MANAGEMENT

The surrounding biotic (i.e., all organisms present in the surrounding habitats) and abiotic (i.e., landscape use, composition and structure) environments of arable fields affect the magnitude of biotic interactions. For example, the crop pollination service provided by native bees depends on the proportion of upland natural habitat surrounding a farm [START_REF] Kremen | The area requirements of an ecosystem service: crop pollination by native bee communities in California[END_REF]. The coexistence of ruderal plants and crops can also increase crop yield by changing the community composition of pollinators [START_REF] Carvalheiro | Natural and within-farmland biodiversity enhances crop productivity[END_REF]. The majority of natural enemies are generalist predators that require non-crop habitats to various extents for reproduction sites, alternative food sources or refuges. A wide range of approaches have been proposed to increase the abundance and diversity of the natural enemy community, including habitat and landscape manipulation and food supplementation [START_REF] Chaplin-Kramer | A metaanalysis of crop pest and natural enemy response to landscape complexity[END_REF]. However, it has proved difficult to assess the contribution of local and landscape factors to pest abundance and to the level of biological control and has led to inconsistent findings [START_REF] Thies | Landscape structure and biological control in agroecosystems[END_REF][START_REF] Bianchi | Sustainable pest regulation in agricultural landscapes: a review on landscape composition, biodiversity and natural pest control[END_REF][START_REF] Chaplin-Kramer | A metaanalysis of crop pest and natural enemy response to landscape complexity[END_REF][START_REF] Veres | Does landscape composition affect pest abundance and their control by natural enemies? A review[END_REF]. Moreover, even though some studies have explored the question of the scale at which biodiversity is important for ecosystem functioning, many gaps remain mainly because many control agents such as arthropods can disperse not only at local scale (i.e., between non-crop habitats and crops) but also at landscape scale (between different crops or different non-crop habitats). A study focusing on two potato pests revealed that predation occurred over different scales for each pest and that its intensity varied with the habitat type: predation on beetle eggs increased in crop and non-crop habitats when field margins were large relative to the potato fields, whereas aphid predation in field margins increased over a small scale with the area of non-crop habitats but did not change in adjacent potato fields [START_REF] Werling | Local and broadscale landscape structure differentially impact predation of two potato pests[END_REF]. Therefore, a better knowledge of these interactions is required to understand changes in predator and parasitoïd communities, non-target effects, effects over different spatial scales (local vs. landscape scales), context-dependent responses and temporal stability of the effects.

AGRICULTURAL OBJECTIVES AND POPULATION DYNAMICS

In agroecosystems, land use distribution has a high turn-over, resulting in strong spatiotemporal variations in food resources for pests. However, a recent meta-analysis revealed that most studies were carried out over a time scale as short as a single sampling season, which is unsuitable for predicting the results of trophic interactions and, therefore, for proposing management strategies [START_REF] Herrera | Land-cover change effects on trophic interactions: current knowledge and future challenges in research and conservation[END_REF]. This is even more likely since a time-lag is often observed in the response of species to predation, competition or parasitism [START_REF] Wangersky | Time lag in prey-predator population models[END_REF] or to changes in land use [START_REF] Jonason | Assessing the effect of the time since transition to organic farming on plants and butterflies[END_REF] and climate [START_REF] Menendez | Species richness changes lag behind climate change[END_REF]. A time-lag is, therefore, to be expected when manipulating biotic interactions through ecological intensification and its duration is likely to depend on a multitude of factors such as species dynamics, the strength of the biotic interaction (generalist vs. specialist) and the landscape context. For instance, a specialist predator will provide efficient control of a prey in the medium term whereas a generalist predator would provide more rapid but less efficient control [START_REF] Hanski | Delayed density-dependence[END_REF][START_REF] Symondson | Can generalist predators be effective biocontrol agents?[END_REF]. The question arises of whether the timing of the response of biological control is compatible with agricultural objectives (e.g., crop yield and farm income) or whether there would be too great an offset in the timing of the control by ecological processes and the change in management system hence limiting potential benefits. Understanding the mechanisms underlying the timing of the response of species after management change is essential for use of ecological intensification by farmers.

ECO-EVOLUTIONARY PROCESSES IN AGROECOSYSTEMS

Organisms are able to adapt to rapid environmental changes and to counteract strong selection pressure. Strong selection and rapid evolutionary changes have been shown in antagonistic interactions such as host-parasite relationships [START_REF] Decaestecker | Host-parasite 'Red Queen' dynamics archived in pond sediment[END_REF][START_REF] Koskella | Advice of the rose: experimental coevolution of a trematode parasite and its snail host[END_REF]. Similarly, rapid natural selection of resistance to pesticides has been shown in various crops and locations (e.g., Rex Consortium, 2013). Consequently, for ecological intensification, the substitution of large-scale chemical selective pressure by a large-scale biotic selective pressure, while avoiding problems of environmental pollution by chemical inputs, could lead to similar problems of pest adaptation. The likelihood that a pest will adapt to its natural enemy and escape the pest control process, which could result in an uncontrolled spread of a pest, is higher if the selective pressure is predictable in time and space [START_REF] Gilligan | Sustainable agriculture and plant diseases: an epidemiological perspective[END_REF]. Including heterogeneity in selection pressures (both in intensity, time and space) could lower the persistence and adaptive potential of pathogens and pests [START_REF] Gilligan | Sustainable agriculture and plant diseases: an epidemiological perspective[END_REF]Rex Consortium, 2013).

The general idea would be to replace one strong selective pressure by several weaker selective pressures carefully distributed in time and space. For example, pathogens could be controlled by weakening the transmission process of a disease by introducing gaps between susceptible and resistant hosts. One way might be to cultivate a mosaic of different crop species, some being susceptible, some being non-susceptible to the disease, or to apply biocides in a mosaic design (see [START_REF] Gilligan | Sustainable agriculture and plant diseases: an epidemiological perspective[END_REF] for depicting the complexity of mechanisms involved in these strategies). An efficient management strategy should thus be drawn up over landscape and crop succession scales since the taxonomic, functional and genetic diversity of the crops are human-managed. Consequently, the nature of the crop (i.e., as a resource or a host) will vary significantly in both time and space in a non-predictable way, reducing the probability of co-adaptation. For management of agroecosystems based on spatio-temporal heterogeneity, the co-evolutionary process of resistance against the enemy and counter-resistance may lead to a mosaic of local adaptations [START_REF] Thompson | The Geographic Mosaic of Coevolution[END_REF], creating de facto gaps in the distribution of pests susceptible and resistant to the enemy (e.g., asynchrony in fluctuations in resistance and susceptibility between demes [START_REF] Thrall | Rapid genetic change underpins antagonistic coevolution in a natural host-pathogen metapopulation[END_REF]. This might be a means of controlling the spread of diseases or pests, provided that the life-cycle and life history traits of both partners are known (e.g., dispersal ability [START_REF] Gilligan | Sustainable agriculture and plant diseases: an epidemiological perspective[END_REF]. For example, resource concentration techniques which are based on the dilution of the host plant in a crop mixture make a parasite or a pest less efficient in locating and colonizing its host plants [START_REF] Ratnadass | Plant species diversity for sustainable management of crop pests and diseases in agroecosystems: a review[END_REF]. The emergence of such geographical mosaics of local adaptation supposes that the spatial heterogeneity is maintained over time. This strategy could not only be applied in agro-forestry and in orchards but also in annual crops to create spatio-temporal mosaics.

The benefits provided by biotic interactions may be outweighed by negative feedback within a network. A recent study of an aphid-transmitted viral pathogen showed that an endophytic fungal mutualist reduced aphid attacks on tall fescue but did not reduce the concentration of the virus, although it did help the plant by alleviating the harmful effects of the viral infection [START_REF] Rua | Fungal endophyte infection and host genetic background jointly modulate host response to an aphid-transmitted viral pathogen[END_REF]. This means that, although plants appear safe by being resistant to a virus, the potential for viral transmission remains the same. This may result in an increase in viral virulence (as shown by [START_REF] Ebert | Evolution -Experimental evolution of parasites[END_REF]: plants with fungi may tolerate higher viral doses, which may be more lethal after their transmission to non-symbiotic plants. Therefore, although it may appear beneficial, protection by a mutualistic symbiont may be detrimental by maintaining the transmission potential of a pathogen. Complex and hidden interactions must be foreseen for successful ecological intensification.

FROM THEORY TO PRACTICE

Ecological intensification can be achieved in different ways. One of the simplest approaches is to replace current intensive practices, for example large monocultures or massive use of pesticides and fertilizers, by intensified biotic interactions without any other major changes in management strategies. This type of intensification based on the use of one or a few species (e.g., natural enemies of a given pest in a context of large monocultures) needs to be thought through in order to avoid pest adaptation in the medium and long term. The most challenging approach is to "intensify" ecological processes in agroecosystems by encouraging natural biodiversity, particularly biocontrol agents, mutualists and symbionts, with the aim of increasing the effectiveness of beneficial biotic interactions. These approaches require redefining farming systems within a holistic approach to agroecosystems to scale up from plot scale to landscape scale and from agricultural timescales to ecological and evolutionary timescales. For ecological intensification to be widely adopted, a set of methods and tools should be available to address the following questions.

(i) What type of biodiversity (e.g., genetic, species, functional) is required to deliver the ecological services (e.g., resources, pests) for the local targets?

The relationship between biodiversity and ecosystem services should be determined for a variety of socio-economic contexts and environmental conditions. Ecosystem functioning is the result of the local functional diversity, the interaction between species and the interaction between species and the local abiotic conditions. It can be translated into ecosystem services for functions that are important for human well-being (Millennium Ecosystem Assessment, 2005;[START_REF] Constanza | The value of the world's ecosystem services and natural capital[END_REF][START_REF] Lamarque | The diversity of the ecosystem services concept and its implications for their assessment and management[END_REF]. This is illustrated by the conceptual framework proposed by Haines-Young and Potschin (2010), which describes the cascade from biodiversity, functions, services, benefits and values to humans. Although the effects of loss of biodiversity on ecosystem functioning have been widely studied over the past two decades, it is still not easy to establish the links between biodiversity and the provision of valuable ecosystem services [START_REF] Cardinale | Biodiversity loss and its impact on humanity[END_REF]. Several functional ecology studies have proposed conceptual frameworks to link biodiversity to ecosystem functions [START_REF] Lavorel | Predicting changes in community composition and ecosystem functioning from plant traits: revisiting the Holy Grail[END_REF][START_REF] Diaz | Incorporating plant functional diversity effects in ecosystem service assessments[END_REF]. For example, abiotic variables and plant traits rather than land use alone were successfully used by [START_REF] Lavorel | Using plant functional traits to understand the landscape-scale distribution of multiple ecosystem services[END_REF] to refine the representation of multiple ecosystem properties and services in a subalpine grassland landscape. Independently, other studies also proposed using diversity indices (e.g., Simpson diversity index, [START_REF] Bateman | Bringing ecosystem services into economic decisionmaking: land use in the united kingdom[END_REF] to quantify the value of biodiversity as a provider of cultural services. Taking account of the true complexity of food webs is one of the difficult challenges in filling the remaining gaps in the understanding of the relationships between biodiversity and ecosystem functioning and services. (ii) How do management practices, local environment conditions and landscape affect biodiversity? Surveys and advanced statistical analyses are required to distinguish the effect of each factor on biodiversity and to examine the interactions between them. Datasets should include data on the biodiversity of interest (abundance of pests, microorganisms, etc), management practices (crop type, tillage, pesticide use, sowing date, etc), local environmental conditions (pH, soil texture, climate, etc) and surrounding biotic (e.g., species richness) and abiotic (landscape composition and structure) environments. [START_REF] Gabriel | Scale matters: the impact of organic farming on biodiversity at different spatial scales[END_REF] recently proposed a novel multi-scale, hierarchical sampling design to account for the effects of farming systems at various scales. However, further study is required to move from the percentage of the variance explained by each group of factors to a full understanding of the ecological processes involved. This could be achieved by translating cropping systems at landscape scale into environmental gradients to allow comparative analysis of very different environmental situations, including situations with very different types of disturbance. Environmental gradients are an excellent way of representing the properties of the environment in order to show how organisms respond to environmental changes [START_REF] Austin | Searching for a model for use in vegetation analysis[END_REF][START_REF] Garnier | Assessing the effects of land-use change on plant traits, communities and ecosystem functioning in grasslands: a standardized methodology and lessons from an application to 11 European sites[END_REF]. This approach was recently proposed at field scale [START_REF] Gaba | Agroecological weed control using a functional approach: a review of cropping systems diversity[END_REF]. This study provides a conceptual framework for characterizing cropping systems on the basis of properties that affect plant performance, i.e., resource and disturbance levels. Cropping system properties have been characterized in three environmental gradients to take account of the system complexity and the relationship between management practices and local conditions, crop types and socio-economic constraints:

(1) the local conditions characterizing the site independently of actual land use, such as climatic conditions, soil and past usage, (2) resource availability, and (3) disturbance levels that depend directly on management practices and types of crop. It should be possible to characterize environmental conditions, cropping systems and landscape composition and structure in the same way. (iii) How, and how rapidly, will organisms adapt to manipulated biotic interactions? This requires detailed knowledge on the co-evolutionary dynamics of the interacting populations. The consequences of the interplay between major processes (biotic interactions such as competition, mutualism and trophic relationships, the capacity of species to adapt to new conditions, etc) should be explored using eco-evolutionary models. Such models have often shown non-trivial effects of interactions between ecological and evolutionary processes on species dynamics and community patterns (e.g., [START_REF] Kearney | Integrating biophysical models and evolutionary theory to predict climatic impacts on species' ranges: the dengue mosquito Aedes aegypti in Australia[END_REF][START_REF] Norberg | Eco-evolutionary responses of biodiversity to climate change[END_REF]. For example, [START_REF] Poisot | Trophic network structure emerges through antagonistic coevolution in temporally varying environments[END_REF] investigated the functioning of trophic enemy-victim networks and suggested that the temporal variation in the productivity of an environment (i.e., the amount of resources available for the prey) is a major factor in co-evolutionary processes.

This is important in agroecosystems, where the resources are maximized to ensure the production of crop biomass. Using this model, together with a geographical structure approach, would be a promising technique for sustainable control of crops' natural enemies. Since the model explicitly takes account of the frequency of resource input, it can be easily transposed for agricultural systems (the crops being the victim and the pest or disease being the enemy) where the frequency of resource input is under human control. (iv) How to ensure the effectiveness and sustainability of ecological intensification (i.e., how long does it take for a biotic interaction to be effective)?

The targeted services must be delivered in a timeframe compatible with farmers' expectations (yield and income) in various pedoclimatic and socioeconomic situations. Markov Decision Process (MDP) or Partially Observed Markov Decision Process (POMDP) have been used to propose new management strategies for conservation [START_REF] Chadès | When to stop managing or surveying cryptic threatened species[END_REF] and agriculture (Tixier et al., 2013b). These mathematical models provide a framework, taking account of the cost-benefit trade-offs, for sequential decision-making when there is a high level of uncertainty. A challenge for ecological intensification in agriculture is to take account of the time-lagged response of populations to biotic interactions in decision making. It is possible to extend this approach to non-Markovian systems taking account of timelagged response [START_REF] Williams | Optimal management of non-Markovian biological populations[END_REF]. A promising approach is to include time-lag in MDP to evaluate whether the application of ecological intensification is compatible with farmers' expectations.

CONCLUSIONS

Although the potential of ecological intensification is easy to understand, there are many constraints related to the complexity of agroecosystems (multiple organisms, diverse temporal, and spatial scales) which prevent this potential being realized. However, current research themes in ecology and evolution are providing useful concepts and tools which improve our ability to gain insights into the functioning of agroecosystems. Such understanding is very important because manipulating biotic interactions is not necessarily gentler than conventional agriculture and may also have undesirable effects such as the introduction of invasive species or loss of biodiversity.

FIGURE 1 |

 1 FIGURE 1 | Picture of a Montagu's harrier (Circus pygargus) hunting a big green grasshopper in a wheat field in the Long Term Ecological Research Network site "Plaine et Val de Sèvre" (Photo: V. Bretagnolle).This shows that agricultural farmlands provide a suitable habitat for species other than the main crop, opening up new possibilities for resource and pest control by using the natural functions of farmland biodiversity.

FIGURE 2 |

 2 FIGURE 2 | Possible means (blue frames) and objectives (red frames) of ecological intensification.

Frontiers

  in Ecology and Evolution | Agroecology and Land Use Systems June 2014 | Volume 2 | Article 29 | 4

ACKNOWLEDGMENTS

Some of the ideas presented were raised during a workshop organized by the INRA, Environment and Agronomy Department. Sabrina Gaba and Laurent Philippot are supported by INRA, Thierry Rigaud by CNRS and François Bretagnolle by University of Burgundy. This work was funded by the ANR project AgrobioSE "ANR-13-AGRO-0001-03." The revised manuscript has been corrected by Tony Tebby.

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Student Member, IEEE Elise C Emond

Alexandre Bousse

Student Member, IEEE Ludovica Brusaferri

Senior Member, IEEE Brian F Hutton

Senior Member, IEEE Kris Thielemans

Improved PET/CT Respiratory Motion Compensation by Incorporating Changes in Lung Density

Keywords: PET/CT, Respiratory motion compensation, lung, density

Positron Emission Tomography/Computed Tomography (PET/CT) lung imaging is highly sensitive to motion. Although several techniques exist to diminish motion artifacts, few account for both tissue displacement and changes in density due to the compression and dilation of the lungs, which cause quantification errors. This paper presents an experimental framework for joint activity image reconstruction and motion estimation in PET/CT, where the PET image and the motion are directly estimated from the raw data. Direct motion estimation methods for motion-compensated PET/CT are preferable as they require a single attenuation map only and result in optimal signal to noise ratio (SNR). Previous implementations, however, failed to address changes in density during respiration. We propose to account for such changes using the Jacobian determinant of the deformation fields. In a feasibility study, we demonstrate on a modified extended cardiac-torso (XCAT) phantom with breathing motion-where the lung density and activity varythat our approach achieved better quantification in the lungs than conventional PET/CT joint activity image reconstruction and motion estimation that does not account for density changes. The proposed method resulted in lower bias and variance in the activity images, reduced mean relative activity error in the lung at the reference gate (-4.84% to -3.22%) and more realistic Jacobian determinant values.

I. INTRODUCTION

P ET lung imaging suffers from patient respiratory motion, affecting the image resolution, localization, and quantification due to the mismatch of the acquired PET data with the attenuation map, typically computed from a single snapshot CT acquisition. A common approach to compensate for respiratory motion, known as motion-compensated image reconstruction (MCIR), is to use image registration to compute deformation fields between different respiratory states (determined using a respiratory signal [START_REF] Nehmeh | Respiratory motion in positron emission tomography/computed tomography: A review[END_REF]), in order to warp the attenuation map and match the PET data, while performing Elise C. Emond was supported by GlaxoSmithKline (BIDS3000030921). Alexandre Bousse was funded by GE Healthcare. Ludovica Brusaferri was funded by Siemens Healthineers and the EPSRC-funded UCL Centre for Doctoral Training in Medical Imaging (EP/L016478/1).

Elise C. Emond, Ludovica Brusaferri, Brian F. Hutton and Kris Thielemans are with the Institute of Nuclear Medicine, University College London, London NW1 2BU, UK (contact: elise.emond.16@ucl.ac.uk).

The Institute of Nuclear Medicine is partially supported by the NIHR University College London Hospitals Biomedical Research Centre.

Alexandre Bousse was with the Institute of Nuclear Medicine, University College London, London NW1 2BU, UK. He is now with LaTIM, INSERM, UMR 1101, Université de Bretagne Occidentale, Brest, France.

iterative reconstruction. This has shown promising results, particularly for improved tumor delineation in lung imaging.

Commonly used MCIR techniques rely on previous determination of the deformation fields from gated images, which are then incorporated into the reconstruction [START_REF] Qiao | A motion-incorporated reconstruction method for gated PET studies[END_REF]. A standard option relies on registering gated PET images reconstructed without attenuation correction [START_REF] Nehmeh | Respiratory motion in positron emission tomography/computed tomography: A review[END_REF] or, more recently, with attenuation correction [START_REF] Lu | Respiratory motion compensation for pet/ct with motion information derived from matched attenuation-corrected gated pet data[END_REF] and accurate time-of-flight (TOF) information [START_REF] Emond | Effect of attenuation mismatches in time of flight PET reconstruction[END_REF], leading to a reduction in localized image artifacts. A second option exploits CINE-CT images to obtain CT-based deformation fields, which should be more robust because of better contrast and higher resolution. However, this is contingent on having corresponding CT and PET gates for good results (therefore requiring amplitude-based respiratory gating). Also, because of the longer CT acquisition, the radiation dose is increased, which cannot be justified in some circumstances.

Instead of deformation fields determined from image registration, the deformation can also be estimated directly from the acquired PET data. Joint reconstruction and motion estimation (JRM) [START_REF] Jacobson | Joint estimation of image and deformation parameters in motion-corrected PET[END_REF]- [START_REF] Bousse | Evaluation of a direct motion estimation/correction method in respiratory-gated PET/MRI with motion-adjusted attenuation[END_REF] is a reconstruction strategy in which, from the entire PET acquisition (as well as a gated respiratory surrogate signal and a single CT image), it is possible to estimate an activity image and the deformation between (respiratory) states. The input single CT image does not need to correspond to one of the gates of the respiratory signal; the computed activity image corresponds to the same respiratory state as the attenuation map (obtained from the CT image). The estimated deformation operator can warp both activity and attenuation images to each gate of the respiratory signal.

The previous methods however usually do not account for the change of density in the lungs, due to respiration. This could be a source of error, especially as both activity concentration and attenuation are affected. In good approximation, the changes can be estimated from local volume changes [START_REF] Cuplov | Issues in quantification of registered respiratory gated PET/CT in the lung[END_REF] (due to overall lung compression and dilation), although the lung mass is not entirely preserved during the respiration. Some fluid exchanges indeed contribute to an additional change of total lung tissue mass [START_REF] El-Chemaly | Lymphatics in lung disease[END_REF], which can be up to 10% over the respiratory cycle [START_REF] Guerrero | Dynamic ventilation imaging from four-dimensional computed tomography[END_REF], where the lowest mass corresponds approximately to end-expiration and the largest to end-inspiration (likely due to the distension of blood vessels following the dilation of the lungs). This would also contribute to the variation of the lung density and tracer concentration during respiration. Accounting for density changes during respiration might be particularly relevant when studying lung diseases such as fibrosis or tumors where the diseased region of the lung is more rigid than other parts of the lung. Therefore, we want to investigate here a method which would comprise a mass-preserving deformation model, which could help with both reducing data inconsistencies and estimating motion more robustly. In this work, the volume changes will be estimated by computing the "Jacobian determinant" (sometimes also referred as "Jacobian") of the deformation field [START_REF] Haber | Numerical methods for volume preserving image registration[END_REF]. Mass-preserving image registration, via incorporation of the Jacobian of the transformation field in the process, has been utilized before [START_REF] Reinhardt | Registration-based estimates of local lung tissue expansion compared to xenon ct measures of specific ventilation[END_REF]- [START_REF] Gorbunova | Mass preserving image registration for lung CT[END_REF], but the idea has not been brought directly into a joint estimation of PET activity image and motion so far.

This paper aims to demonstrate the feasibility of a reconstruction algorithm for respiratory-gated PET data that can take local density changes into account, when only a single CT image is available. We provide an evaluation based on simulated PET data of XCAT phantom images.

In a first section, the theory of mass-preserving registration will be discussed, as well as the modifications necessary to the existing PET JRM to incorporate mass preservation. The method introduced will then be validated using simulated XCAT data of lung acquisition, where density and activity concentration changes are included within the lung.

II. THEORY

This section presents a method built upon [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF]- [START_REF] Bousse | Evaluation of a direct motion estimation/correction method in respiratory-gated PET/MRI with motion-adjusted attenuation[END_REF] and therefore similar notations will be used.

A. Motion-free PET Image Reconstruction

The activity and attenuation maps are respectively denoted λ ∈ C + and µ ∈ C + , where C + = C 0 (R 3 , R + ) is the set of non-negative continuous functions defined on R 3 .

The set of detected prompts g i ∈ N at a detection bin i ∈ 1, n b , where n b is the number of detection bins (either defined from the corresponding detector pair in non TOF case or from the detector pair and the time bin in TOF case), follows a Poisson distribution of expectation ḡi (λ, µ) such that:

g i ∼ Poisson(ḡ i (λ, µ)) , (1)
with ḡi (λ, µ) = τ a i (µ)H i λ + b i , ( 2
)
where τ is the duration of the acquisition, a i (µ) is the attenuation along the line of response L i corresponding to the detection bin i, H i and a i are defined as

H i λ = R 3 λ(r)h i (r) dr (3)
and

a i (µ) = exp - Li µ(r) dr , (4)
where h i : R 3 → R + is the system response function corresponding to the detection bin i, which is compactly supported such that (3) is well-defined, and b i is the number of background events (either scatter or randoms) at the bin i.

The log-likelihood L, omitting the terms independent of λ and µ, is given by:

L(λ, µ) = n b i=1 g i log ḡi (λ, µ) -ḡi (λ, µ) , (5)
where n b denotes the number of detection bins and log(•) the natural logarithm. Assuming the attenuation µ is known (e.g., from CT reconstruction), the activity image λ can be reconstructed by maximizing [START_REF] Jacobson | Joint estimation of image and deformation parameters in motion-corrected PET[END_REF], for example with the expectation-maximization (EM) algorithm [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF].

B. Standard Approach

PET acquisitions are in practice longer than cardiac and respiratory cycles. However, if we assume the acquired PET data can be binned into several n g "gates" such that g = {g } ng =1

where each g ∈ N n b is a vector regrouping the counts corresponding to gate at each bin-during which intra-gate motion can be considered negligible in comparison to PET resolution. The previous model from Section II-A can then be used separately at each gate , [START_REF] Rezaei | Simultaneous reconstruction of the activity image and registration of the CT image in TOF-PET[END_REF] which depends on the activity distribution λ and the attenuation distribution µ at gate through the expectation

∀ ∈ 1, n g , [g ] i = g i, ∼ Poisson(ḡ i, (λ , µ )) ,
ḡi, (λ , µ ) = τ a i (µ )H i λ + b i, , (7)
where τ corresponds to the duration of the binned respiratory gate and b i, is the background term for the -th gate at the detection bin i. The corresponding log-likelihood Λ to be maximized is:

Λ (λ , µ ) = n b i=1 g i, log ḡi, (λ , µ ) -ḡi, (λ , µ ) . ( 8
)
One approach to tackle motion blur is therefore to reconstruct a single image λ 0 from the gated PET data g 0 [START_REF] Bousse | Maximum-likelihood joint image reconstruction and motion estimation with misaligned attenuation in TOF-PET/CT[END_REF]. However this approach suffers from low SNR as only a fraction of the measured counts are used. To address this, each gate can be reconstructed individually, followed by image registration to a single reference gate and averaging, but this approach requires gated attenuation maps (Cine-CT) which is not always available in practice; and increases patient dose. Alternatively, the motion between each of the gates and the reference gate can be estimated (e.g., from CT or magnetic resonance imaging (MRI)) and incorporated in the reconstruction framework; it is referred to as MCIR.

C. Direct Approach

In this section we derive an extended version of the model proposed in [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF]- [START_REF] Bousse | Evaluation of a direct motion estimation/correction method in respiratory-gated PET/MRI with motion-adjusted attenuation[END_REF].

We define the warping operator W ϕ : C + → C + associated to a diffeomorphism ϕ : R 3 → R 3 such that:

W ϕ f = f • ϕ, where f ∈ C + . (9)
In addition, another warping operator can be introduced, using the approximation that the mass is preserved in the lung and both the changes in activity concentration and density during the respiration are directly correlated to the change in volume [START_REF] Cuplov | Issues in quantification of registered respiratory gated PET/CT in the lung[END_REF]. A modified "mass-preserving" warping operator W ϕ can therefore be defined, incorporating into the previous W ϕ the absolute value of the determinant of the Jacobian matrix of

ϕ in R 3 , denoted | det J ϕ |, such that ∀f ∈ C + W ϕ f = | det J ϕ | • W ϕ f , (10)
where the Jacobian matrix,

J ϕ : R 3 → M 3,3 (R) (x, y, z) → ∂ϕ ∂x (r), ∂ϕ ∂y (r), ∂ϕ ∂z (r) , (11)
is differentiable. The definition of this operator is motivated by the consideration of local mass preservation: if we consider a diffeomorphic deformation ϕ which transforms a non-negative continuous distribution

f 1 into f 2 (i.e., f 2 = f 1 • ϕ), then: R 3 f 1 (r) dr = R 3 | det(J ϕ (r)) | f 2 (r) dr , (12)
where det(J ϕ (r)) consequently reflects the volume changes.

At each gate ∈ 1, n g , the corresponding motion is ϕ . The mass-preserving motion-compensated maximumlikelihood is defined by replacing λ and µ by the deformed reference images λ and µ in (8) then summing over the gates:

L(λ, {ϕ } ng =1 , µ) = ng =1 Λ ( W ϕ l λ, W ϕ µ)) . (13)
The addition of mass-preserving conditions within a cost function (not necessarily using |J ϕ |) has been used by various authors, for CT [START_REF] Yin | Mass preserving nonrigid registration of CT lung images using cubic B-spline[END_REF], PET [START_REF] Gigengack | Motion correction in dual gated cardiac PET using mass-preserving image registration[END_REF], [START_REF] Thielemans | Masspreserving image registration using free-form deformation fields[END_REF], single photon emission computed tomography (SPECT) [START_REF] Ue | Nonlinear motion correction of respiratory-gated lung SPECT images[END_REF], [START_REF] Marin | Deformable left-ventricle mesh model for motion-compensated filtering in cardiac gated SPECT[END_REF] and MRI data [START_REF] Chang | A Technique for Accurate Magnetic Resonance Imaging in the Presence of Field Inhomogeneities[END_REF], mostly for image registration. However, mass preservation has never been used, to the best of the authors' knowledge, directly within a motion-compensating PET image reconstruction scheme, using a log-likelihood cost function to estimate the motion between respiratory gates.

D. Discretization

In this work, the discretization is achieved following the scheme proposed in [START_REF] Jacobson | Joint estimation of image and deformation parameters in motion-corrected PET[END_REF], [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF]. Vectors derived from discretization appear in bold. The discretized images corresponding to λ and µ are λ ∈ R nv and µ ∈ R nv respectively, where n v is the number of voxels in the image, and f j = [f ] j and µ j = [µ] j respectively correspond to the activity and attenuation at voxel j ∈ 1, n v .

We use B-spline functions for the deformation field [START_REF] Jacobson | Joint estimation of image and deformation parameters in motion-corrected PET[END_REF], [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF], [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF]. The discrete motion is expressed with an n c control point grid {r n } nc n=1 , r n ∈ R 3 and with a B-spline coefficient vector

α = α X = α X n nc n=1 , α Y = α Y n nc n=1 , α Z = α Z n nc n=1 , where α C n ∈ R is the cubic B-spline coefficient in direction C ∈ {X, Y
, Z} and at control point r n . The corresponding motion ϕ α is defined at each location r ∈ R 3 as on a subgrid of the original image grid comprising n c control points of coordinates {r n } nc n=1 , such that:

ϕ α (r) = r +   nc n=1 α X n B r-rn D nc n=1 α Y n B r-rn D nc n=1 α Z n B r-rn D   ( 14
)
where D is the distance between the control points, r n is the n-th control point and B(r

) = b(x) • b(y) • b(z)
where r = (x, y, z) are the coordinates of a point in the grid and b is the (uniform) cubic B-spline function. We assume here that the control points are equally spaced in x, y and z directions but the implementation could be extended to uneven spacing. The discrete version of the warping operator W ϕ becomes a matrix W α ∈ M nv,nv (R). Additionally, we define the "discrete" determinant ∆ α ∈ R nv by sampling det J ϕα (r) at each voxel.

We redefine the expected number of counts at a gate such as ḡi, (λ, θ, µ) = τ a i ( W α µ) H W α λ i + s i, , where H is the discretized PET system matrix and a i (µ) is the attenuation factor corresponding to the discretized attenuation map µ.

Finally, the discrete version of ( 13) is

L(λ, θ, µ) = ng =1 Λ ( W α λ, W α µ) , (15)
where θ = {α } ng =1 and α is the motion parameter at gate .

E. Joint Image Reconstruction and Motion Estimation

We propose to estimate the activity image λ and the motion parameter ensemble θ by maximizing [START_REF] Yin | Mass preserving nonrigid registration of CT lung images using cubic B-spline[END_REF] with the addition of regularization terms on λ and θ

Φ(λ, θ, µ) L(λ, θ, µ) + βU (λ) + γV (θ) , (16)
where β and γ are the penalty weights associated to the (quadratic) penalty terms U (λ) and V (θ), respectively. Thus, the discussed JRM tries to solve the following optimization problem to estimate an activity image λ and a set of deformation parameters θ:

( λ, θ) ∈ arg max λ, θ Φ(λ, θ, µ) . (17)
The proposed method proceeds by alternating between λ and θ updates, using block sequential regularized expectationmaximization (BSREM) algorithm proposed by Ahn et al. [START_REF] Ahn | Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms[END_REF] for λ-which will not be detailed here-and limitedmemory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) [START_REF]Updating quasi-newton matrices with limited storage[END_REF] for θ-for which we propose a summary in Appendix B. When a mass-preserving model is used, it will be referred as "masspreserving JRM", as compared to "standard JRM".

With the mass-preserving case, the derivation of the gradient of the cost function in the motion update needs to account for the additional dependence in θ due to the addition of Jacobian determinant in the model. The gradient in θ is somehow similar to that of [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF] with the addition of the analytical derivatives of | ∆ α | with respect to the B-spline coefficients α which are obtained in a a similar way as in [START_REF] Yin | Mass preserving nonrigid registration of CT lung images using cubic B-spline[END_REF]. Details on derivation of the cost function are given in Appendix A and B.

F. Image & Motion Estimation Regularization 1) Activity Image regularization:

In standard JRM [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF], a quadratic prior was used to regularize image reconstruction. In this work, U (λ) is the relative difference prior [START_REF] Nuyts | A concave prior penalizing relative differences for maximum-a-posteriori reconstruction in emission tomography[END_REF] and can be expressed as

U (λ) = - nv j=1 v∈Nj ω j,v (λ j -λ v ) 2 (λ j + λ v ) + γ|λ j -λ v | , ( 18
)
where ω j,v in the inverse distance between the center of a voxel j and the center of a voxel v and N j is the neighbourhood of voxel j.

2) Deformation regularization: V (θ) in ( 16) can be expressed as a sum of penalties over all the gates ∈ 1, n g :

V (θ) = ng =1 R(α ) ( 19
)
where R is a penalty on the coefficients {α } ng =1 . In this work, instead of only regularizing the deformation parameters, the regularization is also performed on the Jacobian determinant image. This should prevent the optimization in θ to favor a mass-preserving motion that matches with the noiseas pointed out in [START_REF] Thielemans | Masspreserving image registration using free-form deformation fields[END_REF] for mass-preserving registration-by constraining the Jacobian determinant image to be smooth. The smoothed total variation (STV) penalty [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] was chosen here:

R STV (α ) = - nv j=1 [F | ∆ α |] j 2 2 + ζ 2 , ζ > 0 ( 20
)
where F is the finite forward difference operator used to approximate the gradient of an image and ζ (chosen as 0.3) is used as a smoothing factor, also enabling differentiation in 0. To ensure smooth deformation fields, a small regularization on the B-spline coefficients was also added within V (θ):

R def (α ) = - 1 2 nc n=1 m∈ Nn ω n,m C∈{X,Y,Z} α C ,n -α C ,m 2 , (21)
where ω n,m is the inverse distance between the control points n (i.e., 1, √ 2 or √ 3) and m, N n is the neighborhood of a control point n (within the n c control points grid). The size of the neighborhoods N n was set to 3 × 3 × 3. The values at the image borders were handled by adding padding to the image, using the nearest values.

More details are given in Appendix B.

III. EXPERIMENTS

Five sets of lung XCAT images were generated, corresponding to five different respiratory gates (n g = 5). The corresponding activity and attenuation images are denoted {λ } were varied uniformly in the lung, depending on the total volume change compared to mid-expiration, using XCAT gated lung masks. As a result, compared to the mid-expiration state (gate 3), the changes in density and activity concentration were respectively as follows: gate 1: -11%, gate 2: -6.9%, gate 3: 0%, gate 4: +8.2%, gate 5: +11.5%. It should however be noted that these changes do not correspond to the Jacobian determinant of the true deformation fields from XCAT, therefore the simulations include some model mismatch as would be the case for patient data. For healthy patients, however, we do not expect a uniform change in density and activity concentration throughout the lung, as some regions are expected to ventilate more than others as documented in [START_REF] Guerrero | Dynamic ventilation imaging from four-dimensional computed tomography[END_REF], [START_REF] Cao | Improving intensity-based lung CT registration accuracy utilizing vascular information[END_REF]. The resulting images are shown in Figures 1 and2.

After the activity and activity images were created, the gated projection data were obtained via forward projection of the activity images and separate attenuation correction factors were computed for each gate. Scatters and randoms were added with a level similar to a usual lung acquisition (about 60% of the total number of detected counts). Poisson noise was added to the projection data in the results presented (number of prompts ≈ 3 • 10 7 ), to simulate a standard FDG lung acquisition in our imaging facility. A GE Discovery 710 scanner geometry was used to simulate the data (using GE proprietary projectors in a MATLAB package). For all gates , the background sinogram {s i, } n b i=1 in ( 7) is considered to be known accurately during the entire reconstruction.

The input attenuation map corresponds to the end-expiration state µ 5 , which means that the first gated activity image reconstructed in the algorithm is perfectly aligned with the attenuation map. The number of JRM iterations numJRMIter 0 0.1 was set to 3 (introduced in Appendix B A.).

cm -1 2 3 •10 -2 cm -1
In addition to the mass-preserving JRM presented here, the "standard" JRM (i.e., without mass preservation, as in [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF]) was used to reconstruct the data. A range of regularization values γ for the motion update [START_REF] Gigengack | Motion correction in dual gated cardiac PET using mass-preserving image registration[END_REF] was tested, from weak to strong visual smoothing of the Jacobian determinant images. The image reconstruction parameters were kept identical for all simulations (8 subsets). For each model and penalty weight, 30 different Poisson noise realizations were performed.

The image variance and squared bias in all attenuation images and only in the lung were studied to select the best regularization configuration for each model and to compare the two models. The attenuation image squared bias is given as:

Bias 2 µ = 1 n g ng =1 1 |V | j∈V 1 K K κ=1 µ [κ] j -[µ ] j 2 , ( 22
)
and the attenuation image variance as:

Var µ = 1 n g 1 K -1 ng =1 1 |V | j∈V K κ=1 µ [κ] j -m atn j, 2 , (23)
where

m atn j, = 1 K K κ=1 µ [κ] j , (24)
K = 30, µ [κ] = Ẇ α[κ] µ 1 is the estimated µ image at a
respiratory gate and at the κ-th noise realisation, Ẇ denoting either the mass-preserving JRM warping operator (i.e., Ẇ = W) or standard JRM warping operator (i.e., Ẇ = W). Similarly, we define the image squared bias and variance in the activity images {λ } ng =1 as follows:

Bias 2 λ = 1 n g ng =1 1 |V | j∈V 1 K K κ=1 λ [κ] j -[λ ] j 2 , (25)
and

Var λ = 1 n g 1 K -1 ng =1 1 |V | j∈V K κ=1 λ [κ] j -m act j, 2 , (26)
where κ] is the estimated activity image at a respiratory gate and at the κ-th noise realization and λ [κ] is the MCIR activity image. These measures of variances and squared biases were used to assess the best regularization configuration (using a scatter plot), for each model.

m act j, = 1 K K κ=1 λ [κ] j , (27)
λ [κ] = Ẇ α[κ] λ [
We also define the mean relative error activity image RD ∈ R nv at a gate such that at a voxel i:

[RD ] i = m j (λ ) -[λ ] i [λ ] i . ( 28
)
All measures in the lung were made using the respiratorydependent masks from XCAT, which were eroded using a small kernel.

IV. RESULTS

The scatter plots for the bias-variance study for the attenuation map are given in Figure 3 and for the activity image in Figure 4. The results show that using a mass-preserving model in JRM reduces the image squared biases, especially in the lung region of the activity image. The activity image variances are also decreased.

From the previous bias-variance analysis, we choose to compare the following results using a medium regularization. Reconstructed activity images corresponding to the two different models on θ, are shown alongside the true activity image in Figure 5 (at the reference gate, i.e., end-expiration).

Visually, the differences between the two reconstructed images are not conspicuous. The mean activity ± standard deviation in the lung (using a slightly eroded lung mask given by XCAT at end-expiration), averaged over all noise simulations, was found to be 1559.6 ± 99.7 Bq.mL -1 in mass-preserving JRM and 1533.1 ± 98.7 Bq.mL -1 in standard JRM, where the true activity concentration is of 1612.1 Bq.mL -1 . The mean relative error in the lung mask at end-expiration is -3.22% for mass-preserving JRM and -4.84% for standard JRM. Masspreserving JRM therefore performs better than standard JRM, although the difference is small. When we compute the mean relative error images {RD } ng =1 for both JRM versions, we can see that including mass preservation in the imaging model decreases the error in the lung and outside of the lung over the different gates. The images are shown in Figure 6.

Finally, we want to assess how closely the estimated Jacobian determinant images correspond to the simulated activity concentration and density change ratio images in Figure 7. The values in the liver and in the cardiac region are improved for the mass-preserving model,as they are closer to 1 (this value corresponds to the simulated density and activity concentration change ratio between all gates and end-expiration).

V. DISCUSSION

The results presented in this work show that incorporating a mass-preserving condition within a reconstruction method, such as joint reconstruction and motion estimation (JRM), can reduce quantification errors, especially in the lungs. A direct estimation of the activity image and of the deformation fields between gates can also help with estimating motion more accurately, especially internally in the lung, while needing only one CT acquisition and respiratory-gated PET data.

A bias-variance study over all gates was first done, from a range of penalty weights for the regularization of the deformation parameters. The scatter plots showed that both bias and variance reduce in the gate-dependent activity images (i.e., after using the deformation operator applied to the reconstructed image) when the mass-preserving model is used-both inside and outside the lungs (e.g., in the cardiac region or in the liver). Additionally, the mean relative error images for all different gates were computed. The errors in the lung decreased, especially for the gates which were the most different to the reference end-expiration gate. In particular, for medium regularization on the deformation parameters, the mean relative error in the lung of the reconstructed activity image at end-expiration decreased from -4.84% to -3.22%. The Jacobian determinant images obtained from the deformation fields were also more consistent with the simulated density and activity concentration changes, with values closer to 1 outside of the lungs.

Although we do not show large benefits when the masspreserving method is applied to the XCAT phantom, our approach could have a larger impact for chest acquisitions of patient data. A possible domain of application could be data of patients suffering from diffuse lung diseases, such as idiopathic pulmonary fibrosis or chronic obstructive pulmonary disease; since the lung rigidity varies locally because of the pathology, lung internal motion could be more precisely estimated by considering mass preservation.

A drawback for joint reconstruction in PET is however linked to rather slow convergence and therefore increased computational time (in our settings, using an implementation based on CPU projectors without multi-threading, it took about 11 times longer to process the data with JRM as compared to standard MCIR, i.e., approximately two days), especially when the activity image used to initialize the algorithm has quantification issues (due to misalignment with the µ map). As suggested in [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF], the initialization of the activity image at the start of the alternating strategy is important for both algorithm acceleration and convergence (as the likelihood is not concave, it is necessary to start as close as possible to the true activity image to ensure that it converges). Therefore, it is recommended to initialize the activity image using a gated reconstruction corresponding to the sinogram of the respiratory gate which is the closest to the gate of the input µ map. Additionally, estimating coarse deformation fields between all gates should be beneficial to accelerate the method presented and could also be used for a less noisy activity image initialization using MCIR. To do so, standard methods such as registering non attenuation-corrected gated images could be used [START_REF] Nehmeh | Respiratory motion in positron emission tomography/computed tomography: A review[END_REF]. Another alternative would be the use of a reconstruction of type Maximum-Likelihood reconstruction of Activity and Attenuation (MLAA) [START_REF] Nuyts | Simultaneous maximum a posteriori reconstruction of attenuation and activity distributions from emission sinograms[END_REF], for which the attenuation sinogram can be estimated up to a constant [START_REF] Defrise | Time-of-flight PET data determine the attenuation sinogram up to a constant[END_REF] in TOF PET (within the activity hull). When applied to respiratory-gated data, although MLAA µ maps are expected to be very noisy, recent work on gated data have shown that it could be possible to extract information on motion [START_REF] Hwang | Data-driven respiratory phase-matched attenuation correction without CT[END_REF], when deep learning is used to denoise the images [START_REF] Hwang | Generation of PET Attenuation Map for Whole-Body Time-of-Flight 18 F-FDG PET/MRI Using a Deep Neural Network Trained with Simultaneously Reconstructed Activity and Attenuation Maps[END_REF].

Moreover, a parameterization of the deformation field using cubic B-splines was used in this work. Although B-splines have the advantage of being simple for derivation and computational reasons, other deformation parameterization could be more adapted for lung registration, for example to handle sliding motion against the ribs. Possible approaches could be to extend supervoxel techniques [START_REF] Heinrich | Deformable image registration by combining uncertainty estimates from supervoxel belief propagation[END_REF], [START_REF] Lafitte | Accelerating multi-modal image registration using a supervoxel-based variational framework[END_REF] or to use "fluid material" deformation fields, such as introduced by Christensen et al. [START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF].

Finally, the method presented here could be extended to estimate the values in the µ map jointly with the defor-mation fields and the activity image. Using the additional constraint linked to the modulation of the attenuation and the activity concentration changes with the Jacobian determinant, it may be possible to estimate the attenuation sinograms accurately, determining intrinsically the additive attenuation sinogram [START_REF] Defrise | Time-of-flight PET data determine the attenuation sinogram up to a constant[END_REF] linked to standard MLAA TOF reconstruction.

VI. CONCLUSION

This study shows the feasibility of a mass-preserving joint activity image reconstruction and motion estimation in PET/CT, which could be also extended to PET/MR imaging. Future work will involve validation on patient data.

APPENDIX A: DERIVATION OF THE DETERMINANT OF THE DEFORMATION JACOBIAN MATRIX FOR B-SPLINE REPRESENTATION

The element (k, l) of the 3 × 3 matrix J ϕα (r) is given by:

[J ϕα (r)] k,l = δ k,l + nc n=1 α C k n D ∂B ∂x l r -r n D (29)
where δ k,l is the Kronecker delta for k, l and ∂B ∂x l is the partial derivative of B with respect to x l , C k corresponding to X, Y , Z for k = 1, k = 2, k = 3, respectively, and x l corresponding to x, y or z for l = 1, l = 2, l = 3 respectively. The derivatives of a determinant of a matrix can be obtained via Jacobi's formula, with respect to a given α C k n :

∂ det J ϕα (r) ∂α C k n (α) = tr adj (J ϕα (r)) ∂J ϕα (r) ∂α C k n (α) (30
) where tr(•) is the trace of the matrix and adj(•) is its adjugate. For a given r, the partial derivatives of J ϕα (r) with respect to α C k n are equal to:

∂J ϕα (r) ∂α C k n (α) = 1 D 3 l=1 ∂B ∂x l r -r n D . (31)
By the chain rule, the partial derivative of | det J ϕα (r)| is therefore:

∂ | det J ϕα (r)| ∂α C k n = sgn(det J ϕα (r)) • tr adj (J ϕα (r)) • ∂J ϕα (r) ∂α C k n , (32)
where sgn(•) is the sign function. The non-differentiability of | • | function in 0 is handled by adding a small number to the Jacobian determinant. The analytical derivatives of | det(J ϕα ) | were validated against finite differences in terms of absolute percentage error between the gradient and the finite differences, normalized with respect to the maximum absolute value of the gradient, for random B-spline coefficients and = 0.00001 used in finite differences. The mean and maximum errors were found to be less than 3 • 10 -11 and 4 • 10 -10 respectively.

• Motion Update: maximizing the log-likelihood, from a given activity image, with respect to the B-spline coefficients θ to estimate the deformation between the reference respiratory state (corresponding to the CT image used to compute the input attenuation map) and the PET respiratory gates (binned using a respiratory surrogate signal). The optimization is performed for each respiratory gate separately. • Regularized Image Reconstruction: given the deformation fields estimated during the previous step, MCIR is used to obtain λ, with a modified version of BSREM. Contrary to [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF], the image is always reinitialized to 1 nv before each reconstruction. The algorithm is initialized from a gated regularized image reconstruction at end-expiration (according to the respiratory trace). The first motion estimation is chosen to have more iterations than following ones (here twice as many). A concise pseudo-code is given in Algorithm 1 and some additional details on the motion update and image update are provided in following subsections B and C. A more detailed summary of the implementation used in this work can be found in [START_REF] Bousse | Maximum-likelihood joint image reconstruction and motion estimation with misaligned attenuation in TOF-PET/CT[END_REF].

ALGORITHM 1 Joint Reconstruction of Motion and Activity

Input: µ, gated projection data g θ (0) ← 0 λ (0) ← BSREM(g 1 , µ, β) for r = 1, . . . , numJRMIter do θ (r) ← Motion estimation from (g, µ, θ (r-1) , λ (r-1) , γ) λ (r) ← MC-BSREM(g, µ, θ (r) , β) end for Output: estimated PET image λ, B-spline coefficients θ

B. Motion Update

The penalized log-likelihood Φ(λ, θ, µ) is maximized with respect to θ during the motion update. Let θ (k) be a current estimate of θ at iteration k. θ (k+1) is obtained by performing a maximization along a gradient ascent search direction t (k) :

α (k+1) = θ (k) + δ (k) t (k) (33)
where δ (k) = arg max δ≥0 Φ(λ, θ (k) + δt (k) , µ)

and t (k) is a gradient ascent direction that incorporates the inverse Hessian of Φ, computed with a L-BFGS quasi-Newton line-search algorithm, handling bound constraints [START_REF]Updating quasi-newton matrices with limited storage[END_REF]; although any line-search method could be used. A Fortran implementation [START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF] was used to compute t and δ at each motion estimation iteration. The step length (34) satisfies the Wolfe conditions [START_REF] Moré | Line Search Algorithms with Guaranteed Sufficient Decrease[END_REF], i.e., guarantees sufficient increase of [START_REF] Gigengack | Motion correction in dual gated cardiac PET using mass-preserving image registration[END_REF].

To simplify the notation, dependencies on λ and µ will be omitted in this appendix. For each gate , the expected number of counts at bin i ḡi, ( W α λ, W α µ), are regrouped in a vector ḡ (α ) ∈ R n b . The likelihood L(λ, θ, µ) will be denoted L(θ) and J (•) will represent the Jacobian matrix in α .

Using the chain rule, the gradient of the log-likelihood in α is:

∇ α L(θ) = J (ḡ (α )) ∇ (ḡ (α )) (35)
where ∇ (ḡ (α )) = g /ḡ (α ) -1 R nc . The latter is modified with the incorporation of the masspreserving warping operator W instead of the regular warping operator W used in [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF] (cf Appendix A):

J (ḡ (α )) = -τ diag{H a ( W α µ) W α λ}LJ ( W α µ) + τ H a ( W α µ)J ( W α λ)
where L ∈ M n b ,nv (R) represents the line integral matrix where an element (i, j) is the length of intersection of the line connecting the two detectors corresponding to the bin i with voxel j, H a (µ) = diag {exp (-Lµ)} H and ∀f ∈ {µ, λ}, J ( W α f ) is the Jacobian matrix associated to W α f with respect to α . We have the following formula:

J ( W α f ) = diag{| ∆ α |}J (W α f ) + diag{W α f }J (| ∆ α |)
where J (W α f ) is obtained as in [START_REF] Bousse | Maximum-likelihood joint image reconstruction/motion estimation in attenuation-corrected respiratory gated PET/CT using a single attenuation map[END_REF] and J(| ∆ α |) is obtained by discretizing [START_REF] Hwang | Generation of PET Attenuation Map for Whole-Body Time-of-Flight 18 F-FDG PET/MRI Using a Deep Neural Network Trained with Simultaneously Reconstructed Activity and Attenuation Maps[END_REF].

The gradient ∇Φ, corresponding to the regularized motion estimation step, is equal to ∇L + ∇V , where ∇L is the concatenation of all gate-dependent gradients ∇ α L and ∇V is the gradient of V .

C. Regularized Image Reconstruction

The image update maximizes the motion-dependent penalized log-likelihood (16), using BSREM [START_REF] Ahn | Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms[END_REF], using the relative difference prior introduced in [START_REF] Nuyts | A concave prior penalizing relative differences for maximum-a-posteriori reconstruction in emission tomography[END_REF]. This algorithm was also modified as motion-compensated BSREM (MC-BSREM), using the warping operators introduced in this work. The step size for gated reconstructions (to initialize the joint reconstruction) at a given iteration k was chosen as , in order to diminish the problem of using too large step sizes at the edges of the field of view when parts of the images are only present in some gates. The parameters of the update are the number of subsets numSubs and the penalty weighting factor β; the image update stops either when the median relative change between two iterations is lower than 0.1% inside the XCAT phantom or after 300/numSubs iterations.

- 1 Fig. 1 :

 11 Fig. 1: Modified XCAT activity images at each respiratory gate used for the simulation: original scale (left) and reduced scale (right) for better visualization of the lung activity concentration change, from end-inspiration (top) to end-expiration (bottom); a red horizontal line was plotted to visualize respiratory motion.

Fig. 2 :

 2 Fig. 2: Modified XCAT attenuation images at each respiratory gate used for the simulation: original scale (left) and reduced scale (right) for better visualization of the lung density change, from end-inspiration (top) to end-expiration (right); a red horizontal line was plotted to visualize respiratory motion.

Fig. 3 :

 3 Fig. 3: Tradeoff between variance and the (squared) bias in the µ images, in (a) the entire image and in (b) the lung only, between the two JRM models, using an edge-preserving JRM regularization.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Tradeoff between variance and the (squared) bias in the λ images, in (a) the entire image and in (b) the lung only, between the two JRM models, using an edge-preserving JRM regularization.

Fig. 6 : 5 Fig. 7 :

 657 Fig. 6: Mean relative activity error image at each gate for standard JRM (left) and mass-preserving JRM (right), from end-expiration (top) to end-inspiration (bottom).

1

 1

APPENDIX B: IMPLEMENTATION AND DERIVATION DETAILS

A. Joint Reconstruction Outline

Joint reconstruction of activity image/motion estimation consists of numJRMIter alternations between: