Urszula Ledzewicz

Helmut Maurer
email: maurer@math.uni-muenster.de

Heinz Schättler

OPTIMAL AND SUBOPTIMAL PROTOCOLS FOR A MATHEMATICAL MODEL FOR TUMOR ANTI-ANGIOGENESIS IN COMBINATION WITH CHEMOTHERAPY

Keywords: 2000 Mathematics Subject Classification. Primary: 92C50; Secondary: 49K15, 49M99 Combination therapy, tumor anti-angiogenesis, optimal and suboptimal controls

We consider the problem of minimizing the tumor volume with a priori given amounts of anti-angiogenic and cytotoxic agents. For one underlying mathematical model, optimal and suboptimal solutions are given for four versions of this problem: the case when only anti-angiogenic agents are administered, combination treatment with a cytotoxic agent, and when a standard linear pharmacokinetic equation for the anti-angiogenic agent is added to each of these models. It is shown that the solutions to the more complex models naturally build upon the simplified versions. This gives credence to a modeling approach that starts with the analysis of simplified models and then adds increasingly more complex and medically relevant features. Furthermore, for each of the problem formulations considered here, there exist excellent simple piecewise constant controls with a small number of switchings that virtually replicate the optimal values for the objective.

1. Introduction. Tumor anti-angiogenesis is an indirect cancer treatment approach with the aim to limit the tumor's growth, possibly even shrink the tumor, by depriving it of the vasculature it needs for a steady supply with nutrients. It had already been proposed by J. Folkman in the seventies [START_REF] Folkman | Antiangiogenesis: new concept for therapy of solid tumors[END_REF], but was only enabled with the discovery of inhibitory mechanisms of the tumor in the nineties [START_REF] Folkman | Angiogenesis inhibitors generated by tumors[END_REF][START_REF] Klagsburn | VEGF/VPF: the angiogenesis factor found?[END_REF]. This treatment targets the endothelial cells that form the lining of the newly developing blood vessels and capillaries. These are healthy, genetically stable cell lines and consequently no clonal resistance to angiogenic inhibitors has been observed in experimental cancer [START_REF] Boehm | Antiangiogenic therapy of experimental cancer does not induce acquired drug resistance[END_REF][START_REF] Kerbel | A cancer therapy resistant to resistance[END_REF]. Hence this approach has the advantage of not being susceptible to developing drug resistance [START_REF] Kerbel | Tumor angiogenesis: past, present and near future[END_REF], all too often the limiting factor in chemotherapy. But anti-angiogenic therapy only limits the tumor's support mechanism and thus it seems clear, and any mathematical model and numerous medical studies confirm this, that the tumor will grow back once treatment is halted. Thus it is not efficient as a stand-alone or monotherapy treatment, but it needs to be combined with other mechanisms like chemotherapy or radiotherapy [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF] that kill cancer cells.

Combinations of anti-angiogenic treatment with traditional chemotherapy simultaneously target two compartments, the cancer cells and the vascular cells that support the tumor. Naturally, the question how these treatments should be scheduled arises. For example, chemotherapy needs the vascular system to deliver the cytotoxic agents to the cancer cells while anti-angiogenic treatments precisely target this vasculature. So it might appear reasonable to apply chemotherapy first. On the other hand, tumor angiogenesis is a pathological procedure. With an overexpression of pro-angiogenic factors it creates an intricate, disorganized, and dysfunctional architecture which results in a vascular network with a high fractal dimension. When the tumor is in a stage of uncontrolled growth with the carrying capacity significantly higher than the tumor volume, this vascular network consists of many small capillaries that have a high permeability. As a result, during an allout attack with cytotoxic agents, much of the drug becomes wasted in this elaborate labyrinth of small and leaky vessels and never reaches the primary tumor. Antiangiogenic agents naturally are most effective on these small and weak capillaries and these are the first ones to be destroyed. Consequently, anti-angiogenic treatments have the effect of normalizing the tumor vasculature [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: a new paradigm for combination therapy[END_REF][START_REF] Jain | Vascular normalization as a rationale for combining chemotherapy with antiangiogenic agents[END_REF] and thus enable a much more efficient drug delivery. In the medical literature the word "pruning" has been used in connection with these normalizing capabilities of anti-angiogenic agents. Hence, and there is mounting medical evidence to support this approach, it may be better to start with anti-angiogenic treatment. But these are difficult questions that are far from being answered and the search for an optimal "therapeutic window" when chemotherapy is best applied in relation to anti-angiogenic therapy is of great interest and is currently pursued in medical trials.

In this paper we consider a mathematical model for tumor anti-angiogenesis by Ergun, Camphausen and Wein that has previously been analyzed in combination treatments with radiotherapy [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]. It is a modification of a biologically validated model by Hahnfeldt, Panigrahy, Folkman and Hlatky [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF] who derived a minimally parameterized population based model for tumor growth under antiangiogenic treatment with the primary tumor volume, p, and the carrying capacity of its vasculature, q, as the main variables. The model by Hahnfeldt et al. has been analyzed extensively by d'Onofrio and Gandolfi in [START_REF] Onofrio | Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by[END_REF] and has undergone various modifications and generalizations (e.g., see [START_REF] Onofrio | Rapidly acting antitumoral anti-angiogenic therapies[END_REF][START_REF] Onofrio | A family of models of angiogenesis and antiangiogenesis anti-cancer therapy[END_REF]) including the one considered here. (For modeling related aspects, see the paper [START_REF] Schaettler | Robustness of optimal controls for a class of mathematical models for tumor anti-angiogenesis[END_REF] in this volume.) For the model by Ergun et al., in earlier research [START_REF] Ledzewicz | A synthesis of optimal controls for a model of tumor growth under angiogenic inhibitors[END_REF][START_REF] Ledzewicz | Scheduling of angiogenic inhibitors for Gompertzian and logistic tumor growth models, Discrete and Continuous Dynamical Systems[END_REF], we have given a complete mathematical solution to the monotherapy problem of scheduling an a priori given amount of anti-angiogenic agents in order to maximize the tumor reduction achievable. This is the problem formulation as it was initially posed in [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]; other, related optimal control formulations have been considered by Swierniak (e.g., [START_REF] Swierniak | Modelling combined angiogenic and chemo-therapy[END_REF][START_REF] Swierniak | Direct and indirect control of cancer populations[END_REF]). In [START_REF] Ledzewicz | Bang-bang and singular controls in a mathematical model for combined anti-angiogenic and chemotherapy treatments[END_REF] we have shown that the monotherapy solution indeed becomes a basis upon which the optimal solution for a combination with chemotherapy can be built. More specifically, the same structure of optimal combination therapy solutions is valid as it was verified numerically for a wide range of initial conditions in [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF] for the model by Hahnfeldt et al. [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF]: optimal controls for the anti-angiogenic agent follow the optimal angio-monotherapy and then, at a specific time, chemotherapy becomes active and is given in one full dose session.

In all these papers the simplifying assumption is made to identify the drug dosage with its concentrations. While chemotherapy is administered in the optimal solution in a single maximum dose segment (bang-bang control), the optimal administration of the anti-angiogenic agent follows a time-varying feedback control that depends on the states p and q (singular control). It follows similarly to some of our earlier results for compartmental models for chemotherapy [START_REF] Ledzewicz | The influence of PK/PD on the structure of optimal control in cancer chemotherapy models[END_REF] that optimality of bangbang controls is not affected if a standard linear pharmacokinetic equation is added to the model. But this is different if optimal controls are singular. In this case the so-called order of the singular arc increases [START_REF] Ledzewicz | Singular controls and chattering arcs in optimal control problems arising in biomedicine[END_REF] causing qualitative changes. However, as feedback functions that depend on q, singular controls are not really practical for this problem to begin with and the real question is the following one: what are the effects that the addition of pharmacokinetic equations has on practically realizable and good suboptimal treatment protocols? In this paper, for a fixed set of initial conditions, we give and compare the optimal solutions for four versions of the underlying optimal control problem. We first present the optimal solutions for both the mono-and combination therapy models when drug dosages and their concentrations are identified and then add a standard linear pharmacokinetic equation to model the concentration of the anti-angiogenic agent. Since chemotherapy is given in one full dose session, in order to keep the model simpler, we did not add a second equation for the cytotoxic agent. In fact, this does not significantly change the structure of solutions and here a simplified model that reasonably well reflects the optimal solution is preferred. It will be shown that while the structure of optimal controls becomes quite a bit more complex with this addition, there exist simple realizable suboptimal protocols that provide excellent approximations to the optimal values.

Optimal Control for a Mathematical Model for Combination Therapy.

We formulate a mathematical model for combination of tumor anti-angiogenesis with chemotherapy that is based on a model by Ergun, Camphausen and Wein in [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF]. This model itself is a modification of a biologically validated model by Hahnfeldt et al. from [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF] in which the spatial aspects of the underlying consumption-diffusion processes that stimulate and inhibit angiogenesis are incorporated into a non-spatial 2-compartment model with the primary tumor volume p and its carrying capacity q as variables. Tumor growth is modeled by a Gompertzian growth function,

ṗ = -ξp ln p q , (1)
where ξ denotes a tumor growth parameter. The dynamics proposed in [START_REF] Ergun | Optimal scheduling of radiotherapy and angiogenic inhibitors[END_REF] for the equation modeling the change in the carrying capacity is given by q = bq 2 3 dq

4 3 -µq, (2)
where b (birth) and d (death), respectively, are endogeneous stimulation and inhibition parameters for the endothelial support and the term µq represents natural death terms. The inhibition and stimulation terms, I(q) = dq 4 3 and S(q) = bq 2 3 , are a modification of the corresponding terms, I(p, q) = dp In (2) this dynamics is slowed down by choosing the inhibitory effect of the tumor proportional to the tumor radius and not its surface area as it is done in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF]. Then still replacing p with q in steady state, this results in a significant mathematical simplification of the q-dynamics since the tumor volume p is eliminated from this equation. (For a more detailed description of the underlying modeling, also see the paper [START_REF] Schaettler | Robustness of optimal controls for a class of mathematical models for tumor anti-angiogenesis[END_REF] in this volume).

In combination treatments with chemotherapy two controls u and v are introduced that represent anti-angiogenic and cytotoxic agents, respectively. Antiangiogenic agents typically are biological agents that need to be grown in a lab and still are very expensive and limited. Chemotherapeutic agents, on the other hand, are widely available, but have serious side-effects and thus can only be administered in limited quantities. From a practical point of view it is therefore of importance how given amounts of these agents can be administered to have an "optimal" effect. Mathematically, this can be formulated as an optimal control problem with free terminal time T and isoperimetric constraints that limit the quantities of the agents to be given,

T 0 u(t)dt ≤ A and T 0 v(t)dt ≤ B.
Adding extra variables y and z that monitor the total amounts of agents that have already been administered, and relabeling the total amounts as y max and z max accordingly, mathematically this problem takes the following form:

[C]: for a free terminal time T , minimize the objective J(u) = p(T ) subject to the dynamics

ṗ = -ξp ln p q -ϕpv, p(0) = p 0 , (3) q
= bq 2 3 -dq 4 3 -µq -γqu -ηqv, q(0) = q 0 , (4)
ẏ = u, y(0) = 0, (5)
ż = v, z(0) = 0, (6)
over all Lebesgue measurable functions u : [0, T ] → [0, u max ] and v : [0, T ] → [0, v max ] for which the corresponding trajectory satisfies the end-point constraints y(T ) ≤ y max and z(T ) ≤ z max . The coefficients ϕ, γ and η are non-negative constants that relate the dosages of the respective agents to their effectiveness. The constants u max and v max denote the maximum doses of the anti-angiogenic and cytotoxic agents, respectively, with the total available amounts of each agent denoted by y max and z max . If we formally set z max = 0, i.e., no cytotoxic drugs are available, then we obtain the corresponding anti-angiogenic monotherapy problem and we denote this problem by [M].

Throughout this paper we shall use the following parameter values in our numerical computations: The variables p and q are volumes measured in mm 3 ; ξ = 0.084 per day, b = 5.85 per day, d = 0.00873 per mm 2 per day, These numerical values are based on the data given in [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF]. For illustrative purposes we also chose a small positive value for µ, µ = 0.02 per day. In the formulation [C] dosage and concentration of the agents are identified and we give u and v the meaning of dosages with units of mg of dose per kg per day; the total amounts y max and z max thus are in mg of dose per kg. The conversion factor γ in (4) is taken as γ = 0.15 with units of kg per mg of dose. The maximum dosage for the anti-angiogenic agent is taken to be u max = 15 and we limit the total amount to be y max = 45. These values are chosen for illustrative purposes and are not based on biological data. Regarding the cytotoxic agent we take ϕ = 0.1 kg per mg of dose and we give results for both η = 0 and small positive values for η. For biological reasons the relative effects of the cytotoxic agent on the vasculature should be smaller than on the cancer cells, η ≤ ϕ, and in our computations we only observed the expected quantitative changes. Hence we mostly report numerical results for η = 0. The control limits taken are v max = 20 and z max = 100. We emphasize that our theoretical results about the structure of optimal controls are independent of these parameter values and that they are generally valid. Also, in all our figures we consistently plot p vertically and q horizontally since this better visualizes tumor reductions, the objective in the optimal control problem [C].

Naturally, in this formulation all states need to be non-negative, but this need not be added as an extra constraint. It is not difficult to show that if p 0 and q 0 are arbitrary positive initial conditions, then for any admissible controls u and v, the solution (p, q, y, z) to the corresponding differential equation exists for all times t > 0 and both p and q remain positive (e.g., [START_REF] Ledzewicz | Scheduling of angiogenic inhibitors for Gompertzian and logistic tumor growth models, Discrete and Continuous Dynamical Systems[END_REF][START_REF] Ledzewicz | Bang-bang and singular controls in a mathematical model for combined anti-angiogenic and chemotherapy treatments[END_REF]). It follows from standard results of optimal control (see, e.g., [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF]) that optimal controls for problem [C] exist. But degenerate cases are possible. The reason lies in the fact that, no matter what control is being used, the cancer volume p increases for p < q expϕ ξ v max . Consequently, if the initial condition lies in this region, and if the overall amounts y max and z max simply are too small for the system to reach the region p > q expϕ ξ v max , then indeed the smallest value for p along any solution is given by the initial condition p 0 . Mathematically, in such a case the "optimal" solution simply becomes to do "nothing", i.e., T = 0. Other, less degenerate situations are possible as well in which for similar reasons not all available agents are fully used up. It is not too difficult to analyze the dynamics for the system in these cases as well, but their analysis brings the need to distinguish various subcases. Also, from a practical point of view, these initial data (p 0 , q 0 ) and (y max , z max ) are medically not very realistic. For this reason we make the following definition and henceforth consider an initial condition that is well-posed in this sense. Definition 2.1. We say the initial data (p 0 , q 0 ) and (y max , z max ) are well-posed for problem [C] if the optimal final time T is positive and if all anti-angiogenic and cytotoxic agents have been used up, y(T ) = y max and z(T ) = z max .

A complete solution to the monotherapy problem [M] (i.e., for z max = 0) was given in [START_REF] Ledzewicz | Scheduling of angiogenic inhibitors for Gompertzian and logistic tumor growth models, Discrete and Continuous Dynamical Systems[END_REF]. Optimal controls u * are specific concatenations of the constant controls u = u max (full dose treatment) and u = 0 (no anti-angiogenic agents are being administered) and a specific, time-varying feedback control u sin , a so-called singular control. This control can be calculated explicitly and is given in the Proposition below:

Proposition 1. [20] If the optimal control u * is singular on an open interval I, u * (t) = u sin (t), then u sin (t) = Ψ 3 q * (t) (7)
where

Ψ(w) = b -dw 2 w + 3ξ b + dw 2 b -dw 2 -µ . ( 8
)
This control is only optimal along one specific curve in (p, q)-space, the corresponding singular arc S, which is the graph of a smooth function p sin of q given by p sin (q) = q exp 3 bdq 2 3 µq

1 3 b + dq 2 3 . ( 9
)
The singular control u sin makes this curve invariant. Admissible controls need to satisfy the control constraint 0 ≤ u ≤ u max and this restricts the admissible portion of the singular arc to an interval [q * , q * u ] with the values q * and q * u the unique solutions to the equation

Ψ( 3 √ q) = u max in (0, b d
3 ). Note that u sin (t) is a smooth feedback control that only depends on 3 q * (t). For a typical well-posed initial condition (p 0 , q 0 ) with p 0 ≤ q 0 (i.e., the tumor is growing), optimal monotherapy controls start with an initial segment along which anti-angiogenic agents are given at the maximum dose u max until the singular arc S is reached. At this point optimal controls switch to the singular control u sin and the corresponding optimal trajectory now descends along the singular arc until all anti-angiogenic agents are exhausted. It is here where most of the tumor shrinkage occurs. When the inhibitors have been exhausted, active therapy is over. However, even as no more agents are given, i.e., along a trajectory for the control u = 0, there still is an additional tumor reduction. The reason for this lies in the fact that inhibitors run out in the region p > q where the tumor volume still shrinks, even for u = 0. Hence, because of after effects the minimum tumor volume is only realized when the trajectory corresponding to u = 0 crosses the diagonal p = q. The corresponding time T then is the limit of the horizon considered in the problem formulation [M]. Depending on the actual initial condition (p 0 , q 0 ) and the available amount y max of anti-angiogenic agents, there exist changes in this structure, most notably if the singular control saturates (i.e., reaches its upper limit u max ) before all inhibitors have been used up. But here we concentrate on this most common scenario and refer the reader to [START_REF] Ledzewicz | Anti-angiogenic therapy in cancer treatment as an optimal control problem[END_REF] for a detailed discussion of other aspects of the structure of optimal controls. Fig. 1 shows the synthesis of optimal controlled trajectories for the monotherapy problem [M]. Such a synthesis provides a road map of all optimal trajectories depending on their initial conditions. The center piece is the singular arc S shown as a solid blue curve in the figure. The trajectories corresponding to the constant controls u = 0 (dash-dotted green curves) and u = u max (solid green curves) give the response of the system to these controls and the dotted line marks the diagonal p = q where the maximum tumor reduction is realized. In the diagram several optimal trajectories are indicated and the thick curves in the graph mark one particular such trajectory that follows the concatenation structure u max /u sin /0 described above.

Fig. 2 gives the optimal control and its corresponding trajectory for the initial conditions p 0 = 8, 000 [mm 3 ] and q 0 = 10, 000 [mm 3 ]. In this case the optimal concatenation sequence is u max /u sin /0: the optimal control is given at full dosage u = u max = 15 until the singular curve S is reached at time t 1 = 1.341 days. The administration then follows the time-varying singular control for t 2 = 3.722 days until all anti-angiogenic agents are exhausted after 5.063 days. Due to after effects the maximum tumor reduction is realized along a trajectory for control u = 0 at the optimal terminal time T = 9.378 days when the trajectory reaches the diagonal p = q, u * (t) =    u max for 0 ≤ t < 1.341 u sin for 1.341 ≤ t < 5.063 0 for 5.063 ≤ t < 9.378

. ( 10
)
The theoretically optimal minimum value for these data is given by p [M] = 2242.65. Once chemotherapy is added to the model, this solution to the monotherapy problem actually becomes the basis on which the optimal combination therapy controls are built. Such a property is not generally valid for multi-control problems, but it holds here because of certain relations between higher-order Lie brackets of the drift and control vector fields defining the dynamics (3) and (4). We have the following characterization of the singular control for the anti-angiogenic agent in the presence of chemotherapy.

Proposition 2. [19]

If the optimal anti-angiogenic agent u * is singular on an open interval I, u * (t) = u sin (t), and the optimal cytotoxic agent is given by v * (t) on I, then γu sin (t) + ηv * (t) = Ψ 3 q * (t) [START_REF] Folkman | Angiogenesis inhibitors generated by tumors[END_REF] with Ψ the same function as given above, [START_REF] Onofrio | On optimal delivery of combination therapy for tumors[END_REF]. Furthermore, on the interval I the control v * is bang-bang with at most one switching from v = 0 to v = v max . For v * ≡ 0 the corresponding trajectory must follow the singular arc S of the monotherapy solution; the trajectory is not constrained for v * ≡ v max .

For combination treatments, however, the theoretical analysis of optimal concatenations is more complex and it has not been completed so far. For example, it has not been excluded theoretically yet that the chemotherapeutic agent v could follow a singular regimen. From a biological perspective such a protocol is not welcome since low dosages are known to foster developing drug resistance. In fact, for various mathematical models for cell-cycle specific chemotherapy [START_REF] Ledzewicz | Optimal bang-bang controls for a 2-compartment model in cancer chemotherapy[END_REF][START_REF] Ledzewicz | Analysis of a cell-cycle specific model for cancer chemotherapy[END_REF][START_REF] Swierniak | Optimal control for a class of compartmental models in cancer chemotherapy[END_REF] we have proven that optimal controls are bang-bang. Also, for the model [C] considered here this structure has been supported in many numerical computations. In fact, we have always seen as optimal the following structure: optimal controls for the anti-angiogenic agent u follow the optimal angio-monotherapy and at a specific time τ chemotherapy becomes active and all available cytotoxic agents u are given in one full dose session. For the same initial conditions given above, p 0 = 8, 000 [mm 3 ] and q 0 = 10, 000 [mm 3 ], and small values of η the optimal controls are of the following form

(u * (t), v * (t)) =        (u max , 0) for 0 ≤ t < t 1 (u sin (t), 0) for t 1 ≤ t < t 2 (u sin (t), v max ) for t 2 ≤ t < t 3 (0, v max ) for t 3 ≤ t < T . (12)
Table 1 summarizes some of our numerical results for various values of η. For η = 0 the switching times for the optimal control u * are the same as in the monotherapy case: the control switches from full dose to singular at time t 1 = 1.341 as the singular arc is reached and all anti-angiogenic inhibitors are exhausted at time t 3 = 5.063. Chemotherapy commences earlier at time t 2 = 4.378, but since η = 0 this does not effect the time when anti-angiogenic agents will be exhausted. However, different from the monotherapy case, the minimum tumor reduction no longer is achieved as the trajectory crosses the diagonal. Now we have q(T ) = 2242.65 and the minimum tumor volume is given by p [C,η=0] = 991.14. The reason for this change lies in the fact that the cytotoxic agent gives tumor reductions also in the region p < q which is not possible for monotherapy only. The final value for q(T ) is identical with the final value for the monotherapy problem if η = 0. For positive values of η this naturally no longer is true and the quantitative changes in the optimal values depend on the effectiveness. But overall the structure of the controls does not change. Figs. 3 and4 show the optimal controls and their corresponding trajectories for η = 0.01 and η = 0.02. Note the drop in the dosage of the anti-angiogenic agent u that occurs as the chemotherapy becomes active. It is caused by the presence of the positive ηterm in [START_REF] Folkman | Angiogenesis inhibitors generated by tumors[END_REF] which adjusts the overall effectiveness of the combined anti-angiogenic and cytotoxic effects on the vasculature (see (4)) to follow the optimal relation Ψ as described in Proposition 2.

Mono-and Combination Therapies with a Linear Pharmacokinetic

Model. A commonly made simplification in mathematical models for cancer treatments identifies the agents' dosage with their concentrations and effects. In reality these clearly are different and their relations are studied under the names of pharmacokinetics (PK) and pharmacodynamics (PD). Pharmacokinetic equations model the drugs' concentrations and pharmacodynamics models their effectiveness. The full process is known as drug delivery in the medical literature. The standard and most commonly used model for P K in the pharmaceutical industry is a simple model of exponential growth and decay given by ċ = -kc + u, c(0) = 0, [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: a new paradigm for combination therapy[END_REF] where u denotes the dosage of the agent and c its concentration. The coefficient k is the clearance rate and is related to the half-life of the agents. In this section we show how the optimal controls given above change when this linear pharmacokinetic equation is added to the model. For our numerical computations we take η = 0 in this section.

Clearly, there will be quantitative differences. Equation ( 13) slows down the clearance of the drug and consequently the anti-angiogenic agent will remain effective over a longer time period and thus give lower minimum tumor volumes. But does this change effect the structure of optimal protocols qualitatively? If it doesn't, then simple extensions of our previous results will accommodate this change. If it does, the analysis needs to be modified and possibly done anew. In mathematical models for chemotherapy that we had considered earlier this was not the case [START_REF] Ledzewicz | The influence of PK/PD on the structure of optimal control in cancer chemotherapy models[END_REF][START_REF] Ledzewicz | Optimal controls for a model with pharmacokinetics maximizing bone marrow in cancer chemotherapy[END_REF]. The reason was that optimal controls were bang-bang and this structure is preserved when a linear pharmacokinetic equation ( 13) is added to the model. But optimal solutions for the anti-angiogenic agent are given by singular controls and here this is no longer true. In [START_REF] Ledzewicz | Singular controls and chattering arcs in optimal control problems arising in biomedicine[END_REF] we have shown that while optimality properties of singular arcs are preserved, its so-called intrinsic order increases from 1 to 2. It is well-known in the control literature [START_REF] Bonnard | Singular Trajectories and their Role in Control Theory[END_REF][START_REF] Zelikin | Theory of Chattering Control with Applications to Astronautics, Robotics, Economics and Engineering[END_REF] that an optimal smooth singular control with values in the interior of the control set cannot be concatenated optimally with either of the constant bang controls u = 0 or u = u max . These concatenations are now accomplished by means of chattering arcs, trajectories that correspond to controls that switch infinitely many times between their upper and lower values on a finite interval. This fact is known as the Fuller phenomenon in the optimal control literature. While such a structure may appear odd, it actually is not (the simplest example is a bouncing ball without friction) and it is generic in high dimension [START_REF] Kupka | The ubiquity of Fuller's phenomenon[END_REF]. Thus it needs to be expected and be dealt with. We therefore add the standard linear pharmacokinetic model [START_REF] Jain | Normalizing tumor vasculature with anti-angiogenic therapy: a new paradigm for combination therapy[END_REF] for the anti-angiogenic agent to the mathematical model and, as a simple model for P D, replace the control u in (4) by its concentration c, but otherwise preserve the same formulation. We do not incorporate a second equation for the cytotoxic agent here since we want to demonstrate the effects on one control and it is the anti-angiogenic agent that leads to more interesting phenomena. We thus now consider the following modification of problem [C]:

[CwPK]: for a free final time T , minimize p(T ) subject to ṗ = -ξp ln p q -ϕpv, p(0

) = p 0 , (14) q
= bq 2 3 -dq 4 3 -µq -γqc -ηqv, q(0) = q 0 , (15) ċ
= -kc + u, c(0) = 0, (16) ẏ
= u, y(0) = 0, (17) ż
= v, z(0) = 0, (18)
over all Lebesgue measurable functions u : [0, T ] → [0, u max ] and v : [0, T ] → [0, v max ] for which the corresponding trajectory satisfies the end-point constraints y(T ) ≤ y max and z(T ) ≤ z max .

In this formulation, as before, u and v are dosages measured in mg of dose per kg per day, but now the conversion factor γ in [START_REF] Klagsburn | VEGF/VPF: the angiogenesis factor found?[END_REF] relates the concentration c to the carrying capacity and thus is in units of kg per mg of dose per day.

It is shown in [START_REF] Ledzewicz | Singular controls and chattering arcs in optimal control problems arising in biomedicine[END_REF] that the optimality status of the singular arc is preserved under the addition of such a linear pharmacokinetic model. In fact, the equations that define the singular control and arc in the monotherapy model [M] remain valid, only with a different interpretation. The singular curve S defined by ( 9) is preserved as a vertical surface Ŝ in (p, q, c)-space and the singular arc for the monotherapy problem [MwPK] (i.e., z max = 0) is now given by the intersection of this cylindrical surface Ŝ with the graph of the function c = Ψ 3 √ q defined by ( 8), see Fig. 5. But concatenations with the singular arc are now no longer accomplished by means of the constant controls u = 0 or u = u max , but through chattering controls that switch infinitely often between these values. Fig. 6 gives the graphs of a numerically computed optimal control u for this model and the corresponding concentration c. Again, the same initial conditions p 0 = 8, 000 [mm 3 ] and q 0 = 10, 000 [mm 3 ] were used and in the numerical computations we have taken a value for the half-life k that is supported by experiments with angiostatin, k = 0.38 per day, [START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF]. The optimal control package NUDOCCCS due to Büskens [START_REF] Büskens | Optimierungsmethoden und Sensitivitätsanalyse für optimale Steuerprozesse mit Steuer-und Zustands-Beschränkungen[END_REF] was used to compute a solution of the discretized control problem using nonlinear programming methods with 400 grid points and a seventh order Runge-Kutta integration method. Fig. 6 shows the resulting control and the corresponding concentration. The irregular structure of the numerically computed control on the left is caused by the fact that the theoretically optimal control chatters and has a singular middle segment. As the intervals between the switchings shrink to 0, the numerical values for the controls no longer alternate between their upper and lower values, but already try to connect with the singular control values that lie in the interior of the control interval. The corresponding value of the objective is within the desired stopping criterion for the algorithm. Now the final time is T = 15.548 [days] and the minimum tumor volume is given by p [M,pk] = 564.762 [mm 3 ]. An interesting observation is that this value is significantly smaller than for the model without a pharmacokinetic equation. The reason is that with a clearance rate k the concentration remains high for a much longer time than in the original modeling -we now have T = 15.548 [days] versus T = 9.378 [days] in the model [M] -and the cumulative effects on the carrying capacity q are significant. In some sense the original modeling is that of a fast acting agent that quickly gets cleared. This model clearly favors slow acting agents that have a low clearance rate. Obviously, chattering controls are not of practical interest and the relevant question is what are the effects on the value of the objective. Do there exist good and simple, practically realizable suboptimal protocols? As for the original model [M], this is indeed the case. Fig. 7 gives an example of a suboptimal control that was computed taking the following simple piecewise constant approximation for the optimal control:

ũ(t) =            u max for 0 ≤ t < t 1 0 for t 1 ≤ t < t 2 ū for t 2 ≤ t < t 3 u max for t 3 ≤ t < t 4 0 for t 4 ≤ t ≤ T . ( 19
)
Here both chattering arcs at the beginning and end of the singular arc are approximated by the trajectory of a simple bang-bang control that switches once from u max to 0 and the singular control is approximated by a constant control ū. This particular choice is probably the simplest reasonable approximation to the control structure that the theory predicts as optimal: a chattering control followed by a singular control and one more chattering control. The switching times t i , i = 1, . . . , 4, the final time T , and the value ū of the control are free optimization variables. Using the arc-parametrization method developed in [START_REF] Maurer | Optimization methods for the verification of second order sufficient conditions for bang-bang controls[END_REF] and the code NUDOCCCS [START_REF] Büskens | Optimierungsmethoden und Sensitivitätsanalyse für optimale Steuerprozesse mit Steuer-und Zustands-Beschränkungen[END_REF], we obtain the optimal switching times t 1 = 1.210, t 2 = 3.649, t 3 = 8.294, t 4 = 9.022, the final time T = 15.544 and the constant control ū is given by ū = 3.725. This gives a minimal tumor volume of p[M,pk] = 564.763 that is virtually identical with the optimal minimal tumor volume p [M,pk] = 564.762 for the "chattering control" in Fig. 6. The numerical computations also show that second order sufficient conditions for the underlying optimization problem are satisfied and hence this control is a strict (local) minimum. Fig. 7 shows the optimal suboptimal control on the left and the corresponding concentration on the right. Overall the behavior is very similar as in case of the chattering control, but the system has a much smoother and thus for many aspects preferable response. Since the final time T is much longer, even better approximations are realized in this case than in the monotherapy model [M]. The trajectories for the optimal and sub-optimal controls are indistinguishable. Similar results are valid for the problem [CwPK] with chemotherapy. The structure for the anti-angiogenic inhibitor is the same as for the monotherapy problem with the irregular pattern in the numerically computed control much more pronounced. For these data chemotherapy only becomes active after all inhibitors have been used up. The reason for this is that the P K slows down the clearance of the anti-angiogenic agent to such an extent that the final time T is more than 5 days after all anti-angiogenic inhibitors have been given. Since no P K was included for the chemotherapeutic agent, the time period when it is given at maximum dose, 5 days in this case, is unchanged. Optimization clearly favors to put the chemotherapy last. The final time T = 15.540 is almost identical with the one for the monotherapy problem and the minimum tumor volume is given by p [C,pk] = 249.598. As before, the value of the carrying capacity at the final time, q(T ) = 564.528, is almost identical with the optimal value in the monotherapy problem. In the suboptimal approximation there obviously is no need to approximate the control v that is already bang-bang. Choosing the same structure [START_REF] Ledzewicz | Bang-bang and singular controls in a mathematical model for combined anti-angiogenic and chemotherapy treatments[END_REF] as before, we now obtain the optimal switching times for the anti-angiogenic agent u as t 1 = 1.

Conclusion.

In this paper, for four related mathematical models, we have given optimal solutions for the problem of minimizing the tumor volume with a priori given amounts of anti-angiogenic and cytotoxic agents. For the monotherapy problem [M] this solution is based on a complete analytical solution given earlier [START_REF] Ledzewicz | Scheduling of angiogenic inhibitors for Gompertzian and logistic tumor growth models, Discrete and Continuous Dynamical Systems[END_REF], for the other versions we presented numerically optimal solutions whose structures are in agreement with existing partial theoretical results. Clearly, the results that we presented are for one set of parameters and initial conditions. But we have seen the same structures of the optimal controls when parameters and initial conditions were varied (extreme variations, however, were not considered). For example, in accordance with existing cost limitations the numerical values used here somewhat make the anti-angiogenic agent the limiting quantity, but the structure of the solutions does not change if the total amount y max is increased. Also based on theoretical analysis, there are three important conclusions that come out of our research:

(1) Optimal protocols for the combination therapy models can be built upon the structure of the optimal monotherapy solutions when only anti-angiogenic agents are given. This is not generally true, in fact, quite rare, for arbitrary multi-input nonlinear control systems. It significantly aids in the solution of the combination therapy problem since it allows for an inductive solution procedure.

(2) The addition of a linear pharmacokinetic model changes the qualitative structure of the optimal control u for the anti-angiogenic agent, but only as far as the concatenations with the optimal singular control are concerned. These are now accomplished by means of chattering arcs. As our numerical results show, however, quantitative differences can be substantial.

(3) Although the mathematically optimal solutions are not practical, they provide guidelines on how to find simple, piecewise constant, sub-optimal controls with a very small number of segments that approximate the optimal values arbitrarily closely. These sub-optimal protocols generally can easily be computed as solution to some low-dimensional optimization problem. But it is only the fact that the optimal solution is known that allows to judge their effectiveness.

All our conclusions support a mathematical modeling approach that starts with the analysis of simplified models and then adds increasingly more complex and medically relevant features. If it is possible to build up the solutions step by step -this is the case with the model considered here, but also holds for other mathematical models for tumor anti-angiogenesis -then this approach has the advantage that it clearly brings out the importance of the individual portions (e.g., chemotherapy versus anti-angiogenic treatments). This also is one of the advantages of lowdimensional, minimally parameterized models. Large and complex models that try to fully and accurately reflect the underlying biology necessarily need to rely on simulations and thus only have a limited or local scope. Low-dimensional aggregate models make mathematical analysis possible that can give give more global insights into the underlying mechanisms. For example, the proposed suboptimal protocols are simple and thus realizable. They suggest that few properly timed and chosen constant dosages are an effective means of administering the agents. Naturally, any clinical implementation would be strongly dependent on the accuracy or correctness of the underlying model, but nevertheless our results provide a clear guidance as to what to look for.

Figure 1 .

 1 Figure 1. Optimal synthesis for the model [M]

Figure 2 .

 2 Figure 2. Optimal control (a, left) and corresponding trajectory (b, right) for the monotherapy problem [M] with initial data (p 0 , q 0 ) = (8, 000 mm 3 ; 10, 000 mm 3 )

  η

Figure 3 .Figure 4 .

 34 Figure 3. Optimal controls (left) and corresponding trajectory (right) for the combination therapy model [C] with initial data (p 0 , q 0 ) = (8, 000 mm 3 ; 10, 000 mm 3 ) and η = 0.01. The dosages of the anti-angiogenic and cytotoxic agents are shown on the left in red and blue, respectively.

Figure 5 .

 5 Figure 5. Vertical singular surface Ŝ in (p, q, c)-space and graph of the concentration c = Ψ( 3 √ q)

Figure 6 .Figure 7 .

 67 Figure 6. A numerically computed 'optimal' chattering control u (left) with corresponding concentration c (right) for problem [MwPK]

Figure 8 .

 8 Figure 8. A numerically computed 'optimal' chattering control u (dosages of the anti-angiogenic agent, shown in red on the left) and its corresponding concentration c (on the right) for problem[CwPK]; the control v, the dosage of the cytotoxic agent, is included in blue on the left.

Figure 9 .

 9 Figure 9. A suboptimal bang-bang control u (dosages of the antiangiogenic agent, shown in red on the left) and its corresponding concentration c (on the right) for problem [CwPK]; the control v, the dosage of the cytotoxic agent, is included in blue on the left.

Table 1 .

 1 Optimal switching times for problem [C] for various values of η

		t 1	t 2	t 3	T	p(T )
	0	1.341 4.378 5.063 9.378 991.14
	0.01 1.341 4.604 5.159 9.605 875.15
	0.02 1.341 4.830 5.186 9.830 774.78

  180, t 2 = 3.590, t 3 = 8.286 and t 4 = 8.940. The cytotoxic agent becomes activated at time τ = 10.545 and the final time is given by T = 15.545. The constant approximation for the singular control is, as in the monotherapy problem, ū = 3.725, and the minimum tumor volume is p[C,sub] = 249.598, once again identical with the optimal value.

q and S(p) = bp, chosen in[START_REF] Hahnfeldt | Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy[END_REF]. The main reason for this modification lies in a differential-algebraic nature of the original model with a fast q-dynamics that reaches its steady-state quickly.

Acknowledgement: We would like to thank the referees of this paper for their helpful comments.

This material is based upon research supported by the National Science Foundation under collaborative research grants DMS 0707404/0707410 and DMS 1008209/1008221.


Pierre-Emmanuel Doré
email: pierre-emmanuel.dore@ensta-bretagne.fr

Christophe Osswald
email: christophe.osswald@ensta-bretagne.fr

Arnaud Martin
email: arnaud.martin@univ-rennes1.fr

Continuous belief functions to qualify sensors performances

Keywords: Sensor performances, Continuous belief function, Parametric model, Inference, Fusion

In this paper, we deal with the problem of sensor performance estimation. As we assume that the sensor is described with only few data, we decide to use the theory of belief functions to represent the inherent uncertainty of our information. Hence, we introduce the belief functions framework, especially in the continuous approach. We describe the model of sensor adopted in our study. Knowing the experimental setting, we suggest an approach to model the sources of information describing our sensor. Finally, we combine these sources in order to estimate sensor performances.

Introduction

The theory of belief functions has been introduced by the famous works of A. Dempster about upper and lower probabilities [START_REF] Dempster | Upper and lower probabilities induced by a multi-valued mapping[END_REF][START_REF] Dempster | Upper and lower probabilities generated by a random closed interval[END_REF] and those of G. Shafer [START_REF] Shafer | A mathematical theory of evidence[END_REF] on the theory of evidence. The work of Ph. Smets [START_REF] Smets | Belief functions: the disjunctive rule of combination and the generalized Bayesian theorem[END_REF] contributed to spread it in the scientific community. Recently, thanks to new breakthroughs [START_REF] Strat | Continuous belief functions for evidential reasoning[END_REF][START_REF] Ristic | Belief function theory on the continuous space with an application to model based classification[END_REF][START_REF] Smets | Belief functions on real numbers[END_REF][START_REF] Smets | Kalman filter and joint tracking and classification based on belief functions in the TBM framework[END_REF][START_REF] Caron | Least committed basic belief density induced by a multivariate Gaussian: formulation with applications[END_REF], the application of belief functions on continuous framework has gained some interest. Hence, we apply these results to describe sensor performances.

Wireless sensors networks are more and more used in monitoring applications [START_REF] Meguerdichian | Exposure in wireless ad-hoc sensor networks[END_REF][START_REF] Elfes | Occupancy Grids: A Stochastic Spatial Representation for Active Robot Perception[END_REF]. One crucial issue in this domain is the sensors placement. Indeed, the aim is to place the sensors in order to maximize the chance to detect an intrusion. To fulfil this objective, we have to characterize the performances of a sensor as a detector. To estimate the performance of a sensor, we adopt a parametric approach. As we only have a small amount of measures to define it, we decide to use the belief function framework to take into account the lack of information.

In a first part (section 2), we present the theory of belief function within a continuous frame of discernement. Then, we introduce a parametric model describing sensor performances (section 3). As we have few data, we present some results obtained by A. Dempster and we suggest an approach funded on the maximum necessity and likelihood principles to the estimate the parameters describing the sensor performances (section 4). Hence, considering that the experimental settings provide us three cognitive independant sources of informations (we take measures in three different places), we merge their respective beliefs about the value of the sensor parameters (section 5). Finally, we analyse our results in section 6. Thanks to these operations, we are able to characterize more accurately the sensor performances.

Belief function framework

The theory of belief functions is a tool used to represent the imperfection of a source of information. There are many kinds of imperfection a belief function can describe such as ignorant, vagueness, uncertainty, ... The purpose of this section is to recall some useful parts of this theory.

Belief on real numbers

In [START_REF] Smets | Belief functions on real numbers[END_REF], Ph. Smets describes an approach of belief functions on real number. He suggests to assign mass on the intervals of R R t¡V, Vu. There is a lot of advantage to procede in this way. Indeed, we can easily associate a basic belief density on R to a probability density function on R 2 . However, this framework is quite restrictive and we cannot use it to describe belief function with basic belief assignement on unconnected sets. In a previous work [START_REF] Doré | About using beliefs induced by probabilities[END_REF], we suggest to scan the set of focal elements (the subset of R whose the basic belief assignement is not null), F, using an index function f and a specific index space I.

f I : I ÝÑ F y Þ ÝÑ f I pyq (1)
With this index function, we consider a positive measure µ Ω 1 such as

» I dµ Ω pyq ¤ 1.
Hence, to define a belief function, we have to consider the brace f I , µ Ω ¨. In order to compute belief functions, we need to define for all A in P pΩq (a family of subset of Ω):

F A ty I|f I pyq Au (2)
F A ty I| f I pyq A ¨$ ru (3)
F A ty I|A f I pyqu (4)
Once these subsets are defined, we can compute the following belief functions:

-The belief function:

bel Ω pAq » F A dµ Ω pyq (5)
1 For the sake of simplicity, it happens we considers thatdµ Ω pyq m Ω pf I pyqq.

-The plausibility function:

pl Ω pAq » F XA dµ Ω pyq (6)
-The communality function:

q Ω pAq » F A dµ Ω pyq (7)
In this framework, we define some basic tools. One of them is the conjunctive rule of combination which allow us to merge sources of information. Let

f I 1 , µ Ω 1 ¨and f I 2 , µ Ω 2 ¨two belief functions.
The conjunctive rule of combination [START_REF] Smets | The combination of evidence in the transferable belief model[END_REF] brings the belief function

¡ f I 1 X 2 , µ Ω 1 X 2 ©
[6] such that:

q Ω 1 X 2 pAq q Ω 1 pAq ¤ q Ω 2 pAq (8)
Within this framework, we will study a particular type of belief functions, the consonant ones.

A consonant belief function is a belief function whose the focal sets are nested. This allows us to create a total ordering on F linked to the relation. Hence, we can define an index function f from R , to F such as py ¥ xq ùñ pf pyq f pxqq [START_REF] Smets | Belief functions on real numbers[END_REF]. To generate consonant sets, we can use g, a continuous function from R n to I r0, α max s R . The α-cuts of g are the set such as:

f I cs pαq tx R n |g pxq ¥ αu (9)
We have the property that F cs

A is an element of Borel algebra. Indeed:

F cs A $ r ñ hα inf inf 2 α I|f I cs pαq A $ r @ ñ F cs A |α inf , α max s (10)
Using a similar argument, we can prove that F cs A and F cs A are elements of Borel algebra. Hence, we can define the index function:

f I cs : I r0, α max s ÝÑ 2 f I cs pαq |α I @ α Þ ÝÑ f I cs pαq (11)
If we consider a probability measure µ R n on I, the brace f I cs , µ R n ¨refers to a belief function.

Maximum of necessity

When we work with an "objective" source of information, we can apply the principle of maximum of necessity. This principle comes from the theory of possibility [START_REF] Dubois | The principle of minimum specificity as a basis for evidential reasoning[END_REF][START_REF] Dubois | On possibility/probability transformations[END_REF][START_REF] Mauris | Transformation of bimodal probability distributions into possibility distributions[END_REF]. The theory of possibility is a particular case of the theory of belief functions (focal elements are always nested). In this situation, the plausibility corresponds to possibility distribution and we consider the necessity which

λ(fcs(α1)) dλ(α1) α2 0 dν(α2) α1 λ(fcs(α))=ν([α,αmax]) αmax λ(fcs(α)) dλ(α) α dν(α)
Fig. 1. How to build belief functions from a probability density function Betf .

corresponds to the belief function. The idea is to work with the most informative distribution of possibility (for the necessity ordering) which fulfils the following assumptions. The first one is that the possibility dominates the probability, i.e.

for all A measurable Π pAq ¥ P pAq. The second one is that the ordering must be kept, i.e. P pAq ¥ P pA I q ô Π pAq ¥ Π pA I q.

These conditions can be transposed in the framework of belief functions [START_REF] Mauris | Transformation of bimodal probability distributions into possibility distributions[END_REF]. Finding a belief function which verifies these properties is equivalent to find a nested focal sets family such as for all A belonging to this family, A is the smallest set (for the inclusion ordering) such as P pAq β. This sets family corresponds to the confidence sets in theory of probability. If we have as input a continuous probability density function Betf , the focal set f I cs pα ¨is the α-cuts of Betf . We obtain a belief function defined by f I cs , µ B pRq ¨such as if we adapt the result of [START_REF] Mauris | Transformation of bimodal probability distributions into possibility distributions[END_REF], we obtain [START_REF] Doré | About using beliefs induced by probabilities[END_REF] (cf. figure 1):

pl B pRq pxq 1 ¡ BetP f I cs pBetf pxqq ¨(12)
i.e.2 :

dµ R pαq αdV pαq with V prα, α max sq λ f I cs pαq ¨ [START_REF] Shafer | A mathematical theory of evidence[END_REF] Within this framework, we can build belief functions on real number with complex focal sets. In most of the cases, the information given by a sensor is represented with a probability distribution. However, it happens we do not have enough information to precizely define a distribution a probability. A belief function can be useful to take this phenomenom into account.

Model of sensor

We consider in this paper that a sensor is a detector. Each detector is defined using two caracteristic measures, the probability of true positive (P t ) and the probability of false positive (P f ). The probability of true positive is the probability that we detect an object when there is actually something. The probability of false positive is the probability we decide there is something when there is nothing. In the litterature, the probability of true positive has been defined as a function of the distance between the sensor and the object we want to detect. In our study case, to model the sensor, we make the following assumptions:

1. The sensor is passive. Hence, the probability of false positive is not a function of the distance between the sensor and the object. When the distance between the sensor and the object is growing, the value of P t tends to P f . 2. If the distance between the sensor and the object is smaller than α, the detection is guaranteed. 3. Beyond this point, the probability of true positive decreases such as the inverse of a geometric law of parameter λ. This type of behaviour has been observed on seismic and magnetic sensor or microphone.

Hence, the following equation describes the probability of detection according a distance d between the sensor and the target, the parameters α, λ and P f :

P t pdq 5 1 if d ¤ α α λ ¤ 1 ¡ P f d λ P f otherwise (14)
This is a trade off between the Elfe's model [START_REF] Elfes | Occupancy Grids: A Stochastic Spatial Representation for Active Robot Perception[END_REF] and the geometric model of sensors [START_REF] Meguerdichian | Exposure in wireless ad-hoc sensor networks[END_REF]. The shape of the sensor coverage (the probability of true positive according the distance) described in equation ( 14) is represented in figure 2.

We will estimate the coverage of a sensor, using our parametric model, when we dispose of only a small amount of measures.

Belief functions induced by sampling

The size of the sample set is not large enough to define precizely a distribution of probability. As measure prospecting is an expensive process, we have to find a way to describe imprecise knowlegde on true and false positive probabilities using a small sample set. We will apply our work to the estimation of the sensor performances.

Experimental settings

In order to illustrate our work, we assume we monitor vehicles moving on three roads at different distances of the sensor ( cf. figure 3). In each case, we have As we work on a small amount of data, the estimates of probabilities of true positive are imprecise. In fact, a family of probability distibutions can fit with the data we obtain. The theory of belief functions is an efficient framework to represent a family of probability distributions. Hence, we present some methods to model information induced by a small learning set in the belief functions framework.

The Dempster's approach

The first motive of Dempster [START_REF] Dempster | A Generalization of Bayesian Inference[END_REF] when he developped the theory of upper and lower probabilities was to propose a mathematical framework to model statistical inference. One of his main results in this domain concerns the binomial sampling. When we have a small amount of samples following a binomial law, it is impossible to precizely define the parameter p of this law. Hence we define a set of distributions of probability using belief functions theory. Let P be the interval r0, 1s which describes the parameter p of a binomial law. Let n be the size of the sample set and k be the number of samples belonging to one of the two classes of the binomial sampling. We obtain a conditional belief function m P on the value of p over P knowing n and k [START_REF] Dempster | A Generalization of Bayesian Inference[END_REF] :

m P rn, ks pru, vsq n! pk¡1q!pn¡k¡1q! u k¡1 p1 ¡ vq n¡k¡1 0 k n m P rn, ks pr0, vsq np1 ¡ vq n¡1 k 0 m P rn, ks pru, 1sq nu n¡1 k n (15)
This kind of belief function can model our problem and be used to define our probabilities. Unfortunately, they are difficult to handle and we cannot use on them the max or the min operators. Indeed, generally, these operators do not induce belief functions (cf. example in table 2). Hence, we have to find another way to characterize our information. There are methods using confidence regions or least commitment principle [START_REF] Aregui | Constructing consonant belief functions from sample data using confidence sets of pignistic probabilities[END_REF]. However, in this paper we focus on the maximum of necessity principle as it will sum up the information contained by all the confidence regions in an unique function. Table 2. The max operator does not induce a belief function.

Likelihood and maximum of necessity principles

As we cannot apply max operator on the belief function defined in equation ( 15), we decide to use the maximum likelihood principle [START_REF] Smets | Belief functions: the disjunctive rule of combination and the generalized Bayesian theorem[END_REF] to define a belief function associated to a binomial sampling. This principle can be resumed as following:

Assuming a likelihood function l, the likelihood of a parameter θ knowing that x is true is equal to the likelihood of x knowing that θ is true. lrxs pθq lrθs pxq [START_REF] Smets | Kalman filter and joint tracking and classification based on belief functions in the TBM framework[END_REF] The likelihood function l must satisfy several properties. One of them is that the likelihood of a set of hypothesis is higher than the likelihood of its subsets:

l pΘq ¥ l pθq , θ Θ (17)
Another one is that the likelihood function must be a sub-additive function:

l pθ A θ B q ¤ l pθ A q l pθ B q (18)
Indeed, two disjoint sets of hypothesis can have a likehood equals to 1. Hence, we can use the plausibility as a likelihood function. In our study, according the likelihood principle, we obtain for π subset of P and k in rr0, nss: pl rn, ks pπq pl rn, πs pkq

Hence we have to define pl rn, πs pkq. The distribution of probability P describing a binomial law of parameter p is equal to :

P rn, ps pkq ¢ n k p k p1 ¡ pq n¡k , 0 ¤ k ¤ n (20)
Applying the maximum of necessity principle to this distribution of probability [START_REF] Dubois | On possibility/probability transformations[END_REF], we define a consonant continuous belief function pl rn, ps. This one is entirely described by his contour function computed using algorithm 1. In parametric models, data are highly dependant, hence we prefer not to use a conjunctive rule to merge source of information but a disjunctive one:

pl rn, ks pπq pl rn, πs pkq max pπ ppl rn, ps pkqq

Result: Plausibility of p for a given pair pn, kq. tmp 0; M ¢ n k p k p1 ¡ pq n¡k ; for i r0, ns do P rn, ps piq ¢ n i p i p1 ¡ pq n¡i ; if P rn, ps piq ¤ M then tmp tmp P rn, ps piq; end end pl rn, ps pkq tmp; // maximum of necessity principle. pl rn, ks ppq pl rn, ps pkq; // likelihood principle.

Algorithm 1: How to compute plausibility of p knowing pn, kq.

Combination of information

We have previously defined belief functions to describe the plausibility of P t and P f for given values of d. As we want to define these probabilities on the whole space, we have, using the model of sensor described by equation ( 14), to infer the plausibility of the 3-tuplet pα, k, P f q.

Independant sources of information

Using the sensor model (cf. equation ( 14)), we can associate for each distance d of measure a 3-tuplet pα, k, P f q to a given probability of true positive P t . Moreover, using table 1, we can define pl P d pP t q. Let T be the parameter framework. We obtain thanks to a mass transfert:

pl T d pα, λ, P f q pl P d pP t q (22)
The experiments which lead us to define P t for given distances d are cognitively independent. Hence we can merge these different sources of information using the conjective rule of combination [START_REF] Smets | The combination of evidence in the transferable belief model[END_REF]. As we use consonant belief functions, this operation can be written:

pl T pα, λ, P f q ¹ dt5,10,20u pl T d pα, λ, P f q (23)
Hence, we assume that all the sources of information are relevant and reliable.

Inference using the max

For a given distance d and probability of true positive pP t q, there is, according the sensor model described in equation ( 14), a set D of 3-tuples pα, λ, P f q matching.

As we know the plausibility of each 3-tuples in D, we should be able to compute the plausibility of P t . However, if we use the likelihood principle and set :

pl P d pP t q pl T pDq (24)
we obtain an inextricable equation. Indeed, after combination, the plausibility function described in equation ( 23) is not consonant anymore. Then, we cannot derive the plausibility function from the contour function. Hence, we decide to define a consonant plausibility function such that its contour function fulfil this condition:

pl P d pP t q max pα,λ,P f qD pl T pα, λ, P f q (25)
In this case, this belief function is more specific than the one we would obtain using equation (24). However, as the plausibility of each 3-tuple is built using a singleton P t of P, it can be considered like a good trade-off. The algorithm 2 sums up all the results obtained in section 5 and helps us to compute the plausibity to detect an object at a given distance of the sensor.

Data: cf. table 1

Result: Plausibility function of P t for a given distance d initialization;

for C i pα, λ, P f q r0 : h : 1s ¢ r0 : h : 1s ¢ r0 : h : 1s do for d j t5, 10, 20u do P tpi,djq f unctionpC i , d j q; // cf. equation ( 14)

pl P dj P tpi,djq ¨ f unctionpn dj , k dj , P tpi,djq q; // cf. algorithm 1 pl T dj pC i q pl P dj P tpi,djq ¨;

end pl T pC i q ¹ j pl T dj pC i q; P tpi,dq f unctionpC i , dq; end D 2 C i |abs P tpi,dq ¡ P t ¨¤ h @ ; pl P d pP t q max CiinD pl T pC i q
Algorithm 2: How to compute plausibility.

Results

The results of this work appear in figure 4. The curve d 5m (respectively d 10m and d 20m), represents the plausibility that the probability of true positive of the sensor is equal to P f at the distance d 5m (respectively d 10m and d 20m). To build this curve, we have only used data given by the table 1 and the reasoning described in the section 4. To plot the other curves, we have used the sensor model described in the section 3 and the combination and inference process suggest in section 5.

We remark three things. Firstly, the plausibility function of P t when d 10000 could be linked to the one of P f as P t tends to P f when d increases. Indeed, we have assumed that we use a passive sensor. Secondly, the closer the sensor is of our points of measure, the smaller is the peak of the plausibility function of P t . Hence in this case the estimate of the probability of detection is more accurate. Thirdly, the plausibility functions look like step-wise functions. It is normal as we found the construction of a plausibility function on the binomial law (Bpn, kq) which is discrete.

Conclusion

This study has proposed a way to deal with small amount of data in order to estimate the performance of a sensor. Continuous belief fonctions have been used in two contexts, the first has been to represent an uncertain knowledge about the distribution of probability describing an phenomenon. The second has been to infer information within a parametric model. In this kind of situation, we are often faced inextricable situation and we have to make our approach more simple to find a solution.

This work is a step to define the coverage a sensor networks. Now we have to study the question a the fusion of the information coming from the different sensors in order to have a better estimation of the global network performance. We have also to include some consideration about the reliability of a sensor to improve the study.

Fig. 2 .

 2 Fig. 2. Probability of detection on a grid. The sensor position defines the origin of our coordinates (α 4, λ 0.5, P f 0.15).

Fig. 3 .

 3 Fig. 3. The experimental setting.

plmax

  max ppl1, pl2q 0 0.71 0.59 0.94 0.67 0.93 0.78 1 mmax 0 0.415 0.1723 0.0040 0.1446 0.0923 0.2003 ¡0.0285

Fig. 4 .

 4 Fig. 4. Plausibilities of Pt for different distances.

Table 1 .

 1 in table 1. Generic values.

	Distance sensor/vehicule d number of test n number of detection k
	5 m	20	18
	10 m	20	14
	20 m	20	10

λ refers to the Lebegue's measure.


Paul A Yushkevich

John A Detre

Kathy Z Tang

Angela Hoang2

Dawn Mechanic-Hamilton

María A Fernández-Seara

Marc Korczykowski

Hui Zhang

James C Gee

Hippocampus-Specific fMRI Group Activation Analysis with Continuous M-Reps

Extended Abstract

Continuous M-Reps. In earlier work [START_REF] Yushkevich | Continuous medial representation for anatomical structures[END_REF], we developed a deformable modeling framework called the continuous medial representation (cm-rep), which, in simplest terms, is a continuous analog of the m-rep approach by Pizer et al . [START_REF] Pizer | Deformable m-reps for 3D medical image segmentation[END_REF]. Cm-reps are deformable templates that have a special property that the skeleton (medial axis) of the template preserves its branching topology (number and configuration of its branches) during deformation. Furthermore, the skeleton of the deforming template is modeled explicitly as a continuous manifold and the exact geometric relationship between the skeleton and the template's boundary is captured. The skeleton-boundary relationship makes it possible to define a special shape-based coordinate system on the interior of the template. One of the coordinate axes connects points on the skeleton to the corresponding (nearest) points on the boundary, while the other two axes parameterize the skeleton as well as the boundary. As argued in [START_REF] Yushkevich | Statistical modeling of shape and appearance using the continuous medial representation[END_REF], this coordinate system is a natural way to parameterize the interiors of anatomical structures to which cm-rep templates have been fitted. It allows us to analyze the distribution of appearance features (intensity, fMRI, DTI) within the structures more or less independently of the differences in the shape of the structures. In this paper, we illustrate the application of such analysis to functional MRI data associated with the hippocampus.

fMRI Study. We designed an fMRI study in which hippocampal activation was highly anticipated. Healthy volunteers (n=18) performed a complex scene encoding task in the scanner. The study followed the block design paradigm. In addition to EPI functional imaging, T 1 anatomical scans were acquired. Standard techniques (SPM, [START_REF] Frackowiak | Human brain function[END_REF]) were used to align functional images to the anatomy, to correct for motion artifacts, to smooth the data temporally and spatially and to compute subject-level statistical maps using the general linear model.

Structure-Specific Group Analysis of fMRI. Whole-brain fMRI group analysis methods such as SPM [START_REF] Frackowiak | Human brain function[END_REF] rely on deformable image registration to normalize data between subjects. Registration does not always align the hippocampus accurately [START_REF] Carmichael | Atlas-based hippocampus segmentation in Alzheimer's disease and mild cognitive impairment[END_REF] and such misalignment will, in theory, reduce the statistical significance of intersubject statistical maps. To improve intersubject alignment in the hippocampus, we fit cm-rep templates to manual segmentations of the hippocampus (Fig. 1) and use the cm-rep coordinate system to project subject-level fMRI statistics to a common reference space3 . We then perform random effects analysis to compute a hippocampus-specific groupwise statistical map.

Results: Improved Sensitivity and Specificity. We compared whole-brain and hippocampus-specific maps generated by SPM and our technique, respectively (Fig. 2). The peaks in the hippocampus-specific map had higher raw t-values and lower corrected p-values than the peaks located near the right hippocampus in the whole-brain map (Table 1). This indicates that our approach is more sensitive and that, in general, the accuracy of intersubject normalization does weigh in on the sensitivity of group analysis, despite the influence of many other sources of error. Furthermore, in most studies, peaks and clusters in the whole-brain activation map are assigned anatomical labels by looking up their coordinates in the Talairach space. Due to normalization errors, peaks resulting from consistent hippocampal activation across subjects may be mislabeled as activation in nearby structures. By restricting its analysis to the hippocampus, our approach aims to increase the likelihood that activation labeled as hippocampal is indeed originating in the hippocampus of each subject. However, specificity is confounded by the fact that spatial smoothing and EPI/T 1 alignment are performed for each subject. In part, our future work will focus on developing structure-focused EPI smoothing algorithms. u v boundary skeleton Fig. 1. A cm-rep template fitted to the hippocampus. On the left, the thickness function is plotted in color over the template's skeleton. In the middle, the target segmentation is shown in blue. The shape-based coordinate system is illustrated on the right. Table 1. Peaks near the right hippocampus in the whole-brain group activation map compared to the peaks in the right hippocampus-specific group analysis map. Abbr.: right cerebrum (RC), parahippocampal gyrus (PHG), right hippocampus (RH). The column P hippo indicates the fraction of hippocampus segmentations, warped into the space of the whole-brain template, that overlap each peak's location.

Fig. 2 .

 2 Fig. 2. Examples of sagittal slices through whole-brain (left) and hippocampus-specific (right) group analysis maps. Peaks are identified by circles.

1

  Penn Image Computing and Science Laboratory, Department of Radiology, University of Pennsylvania 2 Center for Functional Neuroimaging, Department of Radiology, University of Pennsylvania

Manual segmentation was performed only for the right hippocampus so our results are currently limited to that structure Mathematical Foundations of Computational Anatomy (MFCA'06)


Brian B Avants

C L Epstein

J C Gee

Geodesic Image Normalization in the Space of Diffeomorphisms

Quantitative, deformable mappings between images are increasingly important as measurement tools in biology and medicine. The theory of diffeomorphisms (smooth differentiable bijections with differentiable inverse) provides a mathematical foundation for the computation and interpretation of these maps. Miller, in particular, has used this theory to develop image normalization techniques that rely on a distance measurement as a regularizer. Here, we use this metric within a general transformation framework that explicitly parameterizes the image to image mapping as a symmetric geodesic path in the space of diffeomorphisms. The key difference between our approach and Miller's is that ours does not choose a single end-point (or template image) from which to measure the map. Rather, we search for the shortest diffeomorphism (smallest deformation) between images by optimizing the transformation with respect to both of its end-points. This algorithm, geodesic normalization, provides solutions that are invariant to which of the input coordinate systems (images) is chosen as a reference. This allows the method to compute metric distances, have truly symmetric performance and to give full space-time solutions that are invertible and diffeomorphic in the discrete domain. Finally, our algorithm guarantees that our solutions and their inverses are consistent to a sub-pixel level.

Introduction

Computational anatomy (CA) uses imaging to make quantitative measurements of the natural world. One may view CA as the science of biological shape and its variation, with roots in the work of Charles Darwin and D'Arcy Thompson [START_REF] Darwin | Origin of Species[END_REF][START_REF] Thompson | On Growth and Form[END_REF]. The wide availability of high resolution in vivo functional and structural imaging has caused a rapid increase in CA's relevance and prominence. Currently, this developing science's primary tools are the topology preserving diffeomorphic transformations. These transformations are used to map an individual image J into the space of a template image, Ī, which serves as a common coordinate system. When a population of images is mapped together by these transformations, each voxel in individual space corresponds smoothly with a single voxel in the reference space. This process creates a continuous spatial map of population information detailing, for example, the relative volume, functional activation or diffusion at a given anatomical position, such as the anterior hippocampus or occipital lobe gray matter. Topology preserving transformation (TPT) are of special interest for this technology as they will introduce neither folds nor tears in this map and will preserve the continuity of curves and surfaces.

Topology preserving normalization permits comparisons to be made across time points in an individual's disease process or to study development patterns across a large population [START_REF] Dale | Cortical surface-based analysis i: Segmentation and surface reconstruction[END_REF][START_REF] Thompson | A surface-based technique for warping 3-dimensional images of the brain[END_REF][START_REF] Csernansky | Hippocampal morphometry in schizophrenia by high dimensional brain mapping[END_REF][START_REF] Studholme | Deformation tensor morphometry of semantic dementia with quantitative validation[END_REF]. Miller, et al. showed that Large Deformation Diffeomorphic Metric Matching (LDDMM) is able to localize hippocampal activation to provide increased statistical significance in functional imaging studies [START_REF] Miller | Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping[END_REF]. Avants et al used a Lagrangian diffeomorphic normalization technique to map and statistically differentiate functionally homologous structures between species [START_REF] Avants | Landmark and intensity-driven lagrangian frame diffeomorphic image registration: Application to structurally and functionally based inter-species comparison[END_REF]. Furthermore, large deformation mappings are better able to separate structural and signal (intensity) differences in population studies, particularly in the presence of atrophy or high shape variation. Diffeomorphic methods may also be extended to normalize vector or tensor images [START_REF] Yan | Large deformation diffeomorphic metric mapping of vector fields[END_REF].

The CA studies cited above were enabled by recently developed, theoretically well-founded methods for studying topology preserving variation. The major advancement in this aspect of CA technology is to base the work in the space of diffeomorphisms. The diffeomorphic space is the broadest smooth, topology preserving space and allows one to very accurately capture both large and small deformation differences in shape. The collection of these transformations forms a mathematical group. Grenander [START_REF] Grenander | General Pattern Theory[END_REF], Mumford [START_REF] Mumford | Pattern theory and vision[END_REF], Miller [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF], Trouve [START_REF] Trouve | Diffeomorphism groups and pattern matching in image analysis[END_REF] and Younes [START_REF] Trouve | On a class of diffeomorphic matching problems in one dimension[END_REF] have studied this group space in the context of computer vision and deformable image transformation and have derived Euler-Lagrange equations for CA [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF].

Our approach to image normalization is based upon Arnol'd's definition of a symmetric, time-parameterized shortest path (geodesic) between two diffeomorphic configurations of a domain [START_REF] Arnold | Sur la gomtrie direntielle des groupes de lie de dimension innie et ses applications l'hydrodynamique des uides parfaits[END_REF]. We argue that this view is fundamental to the theoretical foundation of diffeomorphic normalization, essential for computing true metric distances and desirable for its symmetry properties. Our algorithm will satisfy desirable continuity, anonymity and unanimity conditions [START_REF] Weinberger | On the topological social choice model[END_REF] as well as the metric measurement properties needed for geodesics. This yields a new algorithm, geodesic normalization (GN), that parameterizes the deformation between an image pair with respect both ends of a geodesic path.

Mathematical Background on Diffeomorphisms

We now discuss some basic facts from the theory of diffeomorphisms, the mathematical underpinnings of GN. This section is derived from Arnol'd [START_REF] Arnold | Ordinary Differential Equations[END_REF] and Marsden and Ratiu [START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF]. In this section, we will refer to φ as a geodesic in the space of diffeomorphisms and φ 1 and φ 2 as the components of φ, as described below. We also assume that the images and velocity fields referred to below are sufficiently differentiable and that we are only interested in transformations that visibly change the image. For example, a diffeomorphism of a constant image will Fig. 1. A diffeomorphism is used to map image I into image J (left to center) via coordinate transformation. Composing two diffeomorphisms φ 1 and the inverse of φ 2 enables us to use the maps from K to J and from I to J to make a coordinate transformation from I to K. be reduced to the identity. Image differentiability is required for the derivative computations necessary in the normalization method. Velocity field regularity guarantees the integrability necessary for generating diffeomorphisms. Typically, convolution with a Gaussian ensures image differentiability while a linear operator, L, induces sufficient smoothness on the velocity field. See Dupuis [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF] for a discussion of regularity requirements on images and velocity field and the wellposedness given by diffeomorphic regularization.

Recall that a diffeomorphism is a smooth one-to-one and onto map with a smooth inverse. We always index diffeomorphisms with a spatial coordinate (x, y or z) and, if necessary, a time variable, t. However, we drop the spatial and temporal indexes for brevity where the meaning is clear. We define the diffeomorphic operations that we need as, 1. Coordinate transformation: this operation changes the coordinate system in which an object (image, vector field) is represented. The operation φI transforms the image I(x) into the deformed image Ĩ(y) = I(φ(x, t)) where the intensity I(x) is equivalent to the intensity Ĩ(y). We may also apply this operation to vector (or velocity) fields.

Transformation composition: this operation links diffeomorphisms together,

generating a new map via φ 2 (φ 1 (x, t 1 ), t 2 ).

Understanding these operations are essential in image normalization. We illustrate them in figure 1.

Families of diffeomorphisms can be generated by integrating time-dependent velocity (vector) fields through an ordinary differential equation [START_REF] Arnold | Ordinary Differential Equations[END_REF],

dφ(x, t) dt = V (x, t). (1)
This differential equation defines the change of the map, φ, by the value of a velocity field which is a smooth vector field. Velocities tell us how particles are moving through space: V assigns motion to a specific material point while v defines motion in a fixed coordinate system. The V in equation 1 is a material velocity in the Lagrangian frame. The spatial velocity is computed at the tangent space to the diffeomorphism at time, t, such that v(y, t) = V (x, t) where the change in coordinates is indicated by the use of x in V and y in v. We may also explicitly indicate the change in coordinates by using the map between the two coordinate systems, denoted by φ(x, t) = y, with inverse φ -1 (y, t) = x. Then,

V (x, t) = v(φ(x), t).
Deformable diffeomorphisms commonly used in image registration map domain Ω to itself. The map at the boundary, ∂Ω, may also be defined as the identity, φ(∂Ω) = ∂Ω. This boundary constraint assumes that rigid motion has been factored out of the transformation between images. It also guarantees the transformation is everywhere one-to-one and onto and restricts the solution space to the diffeomorphic subgroup, G 0 .

G 0 is a Frechet Lie group [START_REF] Schmid | Infinite dimensional Lie groups with applications in mathematical physics[END_REF] when G 0 is C ∞ . The length of a diffeomorphic path between elements in this space is similar to the length of a curve, C, connecting two points in Euclidean space, l(C) = 1 0 dC/dt dt, where the Euclidean length of the curve's tangent vector is integrated over its parameterization. Distances in the space of diffeomorphisms are infinite-dimensional analogies of curve length, where the infinitesimal increment in distance is given by a Sobolev norm,

• , operating on the tangent to the diffeomorphism (the spatial velocity) [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF]. A geodesic between ψ 1 and ψ 2 , two elements of G 0 , is defined by taking the infimum over all such paths [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF],

D(φ(0), φ(1)) = inf φ 1 0 v(φ(x, t)) L dt, (2)
φ(0) = ψ 1 and φ(1) = ψ 2 ,
where • L is the Sobolev norm with respect to linear operator, L. Taking the infimum guarantees that we have a geodesic between the elements in G 0 . The length of the geodesic gives a metric distance, does not depend on the origin of its measurement (it has right invariance) and is the basis for GN as well as Miller, Trouve and Younes's work.

A geodesic in the space of diffeomorphisms thus defines the shortest route between two diffeomorphic transformations. Each transformation defines a single, unique configuration of the coordinate system. The length of the path itself is (trivially) symmetric, that is, D(φ(0), φ(1)) = D(φ(1), φ(0) and satisfies metric properties. Furthermore, for all time t ∈ [0, 1], we have φ -1 2 (φ 1 (x, t), 1 -t) = φ 1 (x, 1) = z. Rearranging this equation, we gain intermediate points along the geodesic from φ 2 (z, 1 -t) = φ 2 (φ 1 (x, 1), 1 -t) = φ 1 (x, t). In this way, we see that points along the geodesic are parameterized equivalently from coordinates at either end-point. We will now introduce this coordinate system invariant parameterization into our normalization technology.

Geodesic Image Normalization

The goal of image registration, in general, is to find, for each x in I, the z in J that gives I(x) = J(z) or, alternatively, f (J(z)) where f is an intensityspace transformation. If f is the identity, then the intensity at x in I should be equivalent to the intensity from coordinate z in J. The mapping from x to z may be written φ(x) = z, from image I to image J such that points in I are in one-to-one correspondence with points in J.

When such maps are diffeomorphisms, we include a time parameter, t, that indexes the temporal evolution of φ. We refer to coordinates x in the time zero I domain, z in the time zero J domain and y in a common coordinate system that moves along the curve connecting I(x) and J(z). If this curve is a geodesic, then it is symmetric and will follow the same path whether starting from I or J. This symmetry means our time parameterized maps may be viewed from either endpoint at I or J such that φ 1 (x, t) = y = φ 2 (z, 1 -t). This formulation allows us to deform I and J such that, for any t ∈ [0, 1] I(φ 1 (x, t)) = J(φ 2 (z, 1t)). Such a motion gives a dense map in both space and time and is shown in figure 2 for one point in the image domain. The total mapping between the images is gained through the composition of these two components, φ(x, 1) = φ -1 2 (φ 1 (x, t), 1 -t). GN will exploit this geodesic view to gain symmetry. The algorithm will thus be able to compute the distance between two images, whereas previous algorithms computed asymmetric distance (not symmetric, therefore not a metric distance), due to a biased gradient descent approach that originates in a parameterization of the geodesic with respect to only one endpoint.

Let us now consider the case when we are given two images, I and J, of the same class, known to be (approximately) diffeomorphic. Here, we know neither the path in time nor which image should be considered as the template or reference image. We make the identification of I(0) with I and I(1) with J. We now seek to find the shortest diffeomorphism between these images such that φ 1 (t)I = φ 2 (1 -t)J. We now translate this example into a variational optimization problem. The variational energy for geodesic normalization therefore seeks φ 1 and φ 2 in order to locate the geodesic connecting I and J. Then, φ 1 (x, t)I = φ 2 (z, 1 -t)J, gives the similarity term, |φ 1 (t)I -φ 2 (1 -t)J| 2 . The forward and backward energy is then, using t as a parameter and solving to time t = t,

E GN (I, J) = inf φ1 inf φ2 t t=0 ω{ v 1 2 L + v 2 2 L }dt + Ω |φ 1 ( t)I -φ 2 (1 -t)J| 2 dΩ.
(3) Subject to: each φ i ∈ G 0 the solution of:

dφ i /dt = v i (φ i (t)) with φ i (0) = Id. ( 4
)
Minimization with respect to φ 1 and φ 2 , upholding the arc length constraint, provides the geodesic normalization. Landmarks may also be included in this energy, as in our previous work [START_REF] Avants | Landmark and intensity-driven lagrangian frame diffeomorphic image registration: Application to structurally and functionally based inter-species comparison[END_REF], by dividing the similarity term, as done with the image match terms above. A similar image matching equation appeared in [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF] and [START_REF] Avants | Geodesic estimation for large deformation anatomical shape and intensity averaging[END_REF] as part of a derivation for template generation.

Once this problem is solved, the total symmetric normalization transformation from I to J is φ 1 (x, 1) = φ -1 2 (φ 1 (x, 0.5), 0.5) and from J to I, φ 2 (z, 1) = φ -1 1 (φ 2 (z, 0.5), 0.5). This is distinct from inverse consistent image registration [START_REF] Johnson | Consistent landmark and intensity-based image registration[END_REF] in which a variational term is used to estimate "inverse consistency" and symmetry and invertibility are not guaranteed. The inverse consistency is inherent to our method and is shown for synthetic data in figure 3 and for real data in figure 4. The algorithm is useful for generating shape means as well as symmetric geodesic image interpolation, formulated in [START_REF] Avants | Geodesic image interpolation: Parameterizing and interpolating spatiotemporal images[END_REF].

We will show below that the maps computed by this algorithm satisfy Eckmann's continuity, symmetry (anonymity) and unanimity conditions. Eckmann and Weinberger discuss the existence of such maps [START_REF] Weinberger | On the topological social choice model[END_REF][START_REF] Eckmann | Social choice and topology[END_REF] and note the connection with the generalized mean. Denote a symmetric map, φ, connecting I and J. Eckmann's properties, adapted for image normalization, are then 1. Continuity: the map should vary continuously with the inputs I, J. 2. Symmetry / Anonymity: φ does not depend on permutations of I, J. This is verified if A(I, J) = φ then A(J, I) = φ -1 . 3. Unanimity [START_REF] Weinberger | On the topological social choice model[END_REF]: the map should output φ = Id if I = J.

A map that violates anonymity (2 above) is labeled asymmetric. We now argue that GN satisfies the three conditions given above.

Theorem 1. The solutions, φ 1 (x, t) and φ 2 (y, s), found by GN satisfy the three generalized mean axioms above.

Continuity. Continuity was shown in Dupuis's proof of well-posedness for the diffeomorphic variational image matching problem [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF]. As our problem is, in a global sense, identical, continuity is inherited. Unanimity. The unanimity condition is also satisfied as all velocities will be 0 if all images are identical. Anonymity. We prove anonymity by simply checking that permuting the labels in the variational energy produces an identical optimization problem, which visual inspection confirms. Furthermore, the Euler-Lagrange equations depend symmetrically on gradients of both I and J.

. Geodesic Normalization Implementation: Algorithm 1 states the locally optimal GN algorithm without landmarks, as described in previous sections. A similar approach can be used for landmark matching. The input to the algorithm is a pair of images and a user set number of time interpolation points, {t i }, where t 0 = 0 and t n = 0.5 and n is even. The output is φ 1 , φ 2 defined on Ω × [0, 0.5] and their inverses. These define the full φ i mappings and their inverses at (deformation) time 1.

Algorithm 1 : Geodesic Normalization (GN)

The algorithm notation is the same as in the body of the paper. We compute φ -1

1 and φ -1 2 with the inversion method used in our Lagrangian Push Forward (LPF) method [START_REF] Avants | Landmark and intensity-driven lagrangian frame diffeomorphic image registration: Application to structurally and functionally based inter-species comparison[END_REF].

1:

∀t i Initialize φ 1 (t i ) = φ 2 (t i ) = Id. 2: while ∆E GN > 1 do 3:
Set E old = E new .

4:

Compute φ 1 (x, 0.5)I and φ 2 (x, 0.5)J.

5:

Compute velocities from the symmetric Euler-Lagrange equations of 3.

6:

Set v 1 = v 2 = min( v 1 , v 2 ).

7:

Update φ 1 (0.5) and φ 2 (0.5) by gradient descent such that, for i = 1, 2, φ i (x, 0.5) = φ i (x, 0.5) + λV i (x, 0.5) where V (x, 0.5) = v(φ i (x, 0.5)) and λ is a gradient step length.

8:

Starting at time 0.5 and going backwards toward 0, for all time points, t i where 0 < t i < 0.5, update φ 1 , φ 2 by gradient descent on the length of φ 1 and φ 2 . Note that φ i (0) = Id and φ i (0.5) does not change in this step.

9:

Use the inversion method (below) to find φ -1 1 and φ -1 2 over all time.

10:

Compute E new from equation 3.

11:

∆E GN = E old -E new . 12: end while Symmetry (guaranteed sub-pixel invertibility and algorithmic independence to input permutations) is built into GN and allows us to symmetrically match images to the degree that discrete diffeomorphisms are invertible. An example of a symmetric image registration when large deformations are present is shown in figure 3. Here, we use GN to find spatial correspondences by deforming a half C to a full C. We also illustrate GN for normalizing severely atrophied brains, some of which also suffer from the presence of lesions, in figure 4.

Discussion

A common problem with image normalization algorithms, in general, is asymmetry. Registering image I to J may not produce the same correspondences as registering J to I. This problem has been addressed by methods that compute mappings in both directions. Thirion developed such an approach for Demons image registration [START_REF] Thirion | Image matching as a diffusion process: an analogy with maxwell's demons[END_REF]. Christensen's inverse consistent image registration (ICIR) [START_REF] Christensen | Consistent image registration[END_REF] uses a similar idea. Both algorithms rely upon estimating a measure of "consistency", defined as the difference between the mapping from I to J and the mapping from J to I. Define x in image I, its displacement u(x) and their sum as y = x + u(x). Similarly, define z in image J and its displacement w(z). Consistency at x is C(x) = y + w(y) -x = u(x) + w(y) . Both algorithms attempt to minimize C(x) and, similarly, C(z). However, neither method guarantees the inverse's existence nor gives a well-defined numerical method for its computation. Geodesic image normalization subsumes the above approaches by formulating the normalization process in space and time with a diffeomorphic parameterization. We do not have to "check" the consistency as (via Theorem 1) the consistency is guaranteed by the sub-pixel invertibility of φ 1 and φ 2 and the fact that these transformations may be composed together. Furthermore, our methods explicitly optimize these transformations in the large deformation space.

Conclusion

We explicitly optimize the length of φ, parameterized symmetrically by its components, φ 1 , φ 2 . Optimizing the length of a diffeomorphism is an alternative to finding a geodesic by solving the Euler equations [START_REF] Arnold | Topological methods in hydrodynamics[END_REF]. Our method is shown to satisfy the axioms of continuity, symmetry and unanimity. Satisfying these properties eliminates problems of algorithmic asymmetry. Finally, GN gives a robust estimate to geodesic distance that is invariant to which image, of a pair, is selected as template. Future work will focus on empirically demonstrating improved performance due to symmetrically parameterizing the large deformation image normalization problem. Fig. 4. We volumetrically map an average image template to a set of lesioned and/or atrophied brains using GN. The left column shows a slice of the original subject image. The second column shows the template mapped to the subject space. The center column shows the subject mapped to template space. The second to last column shows the grid deformation from template to subject. The final column shows the original template slice. Asymmetric methods do not perform as well for normalizing this dataset, particularly in the presence of difficult and unpredictable lesions.

Fig. 2 .

 2 Fig. 2. An illustration of the geodesic path taken, in time, by the particle at position x in domain Ω. The path is geodesic if the diffeomorphism associated with the domain minimizes the distance metric in equation 2. The geodesic, φ is the whole path. Its points are traversed via φ 1 and φ 2 .

Fig. 3 .

 3 Fig. 3. The 1 2 C (a) is registered to the full C (d). The grid of the deformation from (a) to (d) is shown in figure (b). The Jacobian of the transformation from C to 1 2 C is in (c) while its histogram is in (f). The 1 2 C to C result is in (e).