Jirka Maršík
email: jirka.marsik@loria.fr

Maxime Amblard
email: maxime.amblard@loria.fr

Introducing a Calculus of Effects and Handlers for Natural Language Semantics

Keywords: compositionality, side effects, monads, handlers, deixis, conventional implicature

. We prove that our calculus exhibits some of the key formal properties of the lambda calculus and we use it to construct a modular semantics for a small fragment that involves multiple distinct semantic phenomena.

Introduction

The prevailing methodology of formal semantics is compositionality in the sense of Frege: denotations of complex phrases are functions of the denotations of their immediate constituents. However, several phenomena have been identified that challenge this notion of compositionality. Examples include anaphora, presupposition, quantification, deixis and conventional implicature. In all of these examples, simple models of denotation (i.e. noun phrases are individuals, sentences are truth-values) run into complications as the denotations can depend on external values (anaphora, deixis) or on something which is not an immediate constituent (presupposition, quantification, conventional implicature).

Among the solutions to these challenges, we find (at least) two types of solutions. First, we have those that relax the condition of compositionality. Notably, the denotation of a complex phrase is no longer a function per se of the denotations of its immediate subconstituents. Rather, it is some other formally defined process. 1 Examples of this approach include:

the incremental algorithm used to build discourse representation structures in DRT, as presented in [START_REF] Kamp | From discourse to logic[END_REF] the λµ calculus, used in [START_REF] De Groote | Type raising, continuations, and classical logic[END_REF] to analyze quantification, since, due to the lack of confluence, function terms do not denote functions over simple denotations the use of exceptions and exception handlers in [START_REF] Lebedeva | Expression de la dynamique du discours à l'aide de continuations[END_REF] to model presuppositions in an otherwise compositional framework the parsetree interpretation step in the logic of conventional implicatures of [START_REF] Potts | The logic of conventional implicatures[END_REF] that builds the denotation of a sentence by extracting implicatures from the denotations of its subparts (including the non-immediate ones)

The other approach is to enrich the denotations so that they are parameterized by the external information they need to obtain and contain whatever internal information they need to provide to their superconstituents. Here are some examples of this style:

any kind of semantic indices (e.g. the speaker and addressee for deixis, the current world for modality), since they amount to saying that a phrase denotes an indexed set of simpler meanings the continuized semantics for quantification [START_REF] Barker | Continuations and the nature of quantification[END_REF] in which denotations are functions of their own continuations • and more generally, any semantics using type raising or generalized quantifiers for noun phrase denotations the dynamic denotations of [START_REF] De Groote | Towards a montagovian account of dynamics[END_REF] that are functions of the common ground and their continuation compositional event semantics, such as the one in [START_REF] Qian | Event in compositional dynamic semantics[END_REF], that shift the denotations of sentences from truth-values to predicates on events

We want to find a common language in which we could express the above techniques. Our inspiration comes from computer science. There, a concept known as monad has been used:

in denotational semantics to give the domain of interpretation for programming languages that involve side effects [START_REF] Moggi | Notions of computation and monads[END_REF]. in functional programming to emulate programming with side effects via term-level encodings of effectful programs [START_REF] Wadler | The essence of functional programming[END_REF].

These two principal applications of monads align with the two approaches we have seen above. The one where we change our calculus so it no longer defines pure functions (e.g. is non-deterministic, stateful or throws exceptions) and the one where we use a pure calculus to manipulate terms (denotations) that encode some interaction (e.g. dynamicity, continuations or event predication). Monad is a term from category-theory. Its meaning is relative to a category. For us, this will always be the category whose objects are types and whose arrows are functions between different types. A monad is formed by a functor and a pair of natural transformations that satisfy certain laws. In our case, this means that a monad is some type constructor (the functor part) and some combinators (the natural transformations) that follow some basic laws. To give an example of this, we can think of the functor T (α) = (α → o) → o together with combinators such as the type raising η(x) = λP. P x as a monad of quantification.

The relationship between side effects in functional programming and computational semantics has been developed in several works [START_REF] Shan | Linguistic side effects[END_REF][START_REF] Van Eijck | Computational semantics with functional programming[END_REF],2 stretching as far back as 1977 [START_REF] Hobbs | Making computational sense of montague's intensional logic[END_REF]. The usefulness of monads in particular has been discovered by Shan in 2002 [START_REF] Shan | Monads for natural language semantics[END_REF]. Since then, the problem that remained was how to compose several different monads in a single solution. Charlow used the popular method of monad morphisms3 to combine several monads in his dissertation [START_REF] Charlow | On the semantics of exceptional scope[END_REF]. Giorgolo and Asudeh have used distributive laws to combine monads [START_REF] Giorgolo | Natural language semantics with enriched meanings[END_REF], while Kiselyov has eschewed monads altogether in favor of applicative functors which enjoy easy composability [START_REF] Kiselyov | Applicative abstract categorial grammars[END_REF].

Our approach follows the recent trend in adopting effects and handlers to combine side effects [START_REF] Bauer | Programming with algebraic effects and handlers[END_REF][START_REF] Kammar | Handlers in action[END_REF] and to encode effectful programs in pure functional programming languages [START_REF] Kiselyov | Extensible effects: an alternative to monad transformers[END_REF][START_REF] Brady | Programming and reasoning with algebraic effects and dependent types[END_REF].

The idea is that we can represent each of the relevant monads using an algebra. We can then combine the signatures of the algebras by taking a disjoint union. The free algebra of the resulting signature will serve as a universal representation format for the set of all terms built from any of the source algebras and closed under substitution. Then, we will build modular interpreters that will give meanings to the operators of the algebras in terms of individuals, truth-values and functions.

In Sect. 2, we will introduce a formal calculus for working with the algebraic terms that we will use in our linguistic denotations. In Sect. 3, we will incrementally build up a fragment involving several of the linguistic phenomena and see the calculus in action. Before we conclude in Sect. 5, we will also discuss some of the formal properties of the calculus in Sect. 4.

Definition of the Calculus

Our calculus is an extension of the simply-typed lambda calculus (STLC). We add terms of a free algebra into our language and a notation for writing handlers, composable interpreters of these terms. An operator of the free algebra corresponds to a particular interaction that a piece of natural language can have with its context (e.g. a deictic expression might request the speaker's identity using some operator speaker in order to find its denotation). A handler gives an interpretation to every occurrence of an operator within a term (e.g. direct speech introduces a handler for the operator speaker that essentially rebinds the current speaker to some other entity).

Having sketched the general idea behind our calculus, we will now turn our attention to the specifics. We start by defining the syntactic constructions used to build the terms of our language.

Terms

First off, let X be a set of variables, Σ a typed signature and E a set of operation symbols. In the definition below, we will let M , N . . . range over terms, x, y, z. . . range over variables from X , c, d. . . range over the names of constants from Σ and op, op i . . . range over the operation symbols in E.

The terms of our language are composed of the following:

M, N ::= λx. M [abstraction] | M N [application] | x [variable] | c [constant] | op M p (λx. M c ) [operation] | η M [injection] | op 1 : M 1 , . . . , op n : M n , η: M η N [handler] | -• M [extraction] | C M [exchange]
The first four constructions -abstraction, application, variables and constants -come directly from STLC with constants.

The next four deal with the algebraic expressions used to encode computations. Let us sketch the behaviors of these four kinds of expressions.

The operation (op) and injection (η) expressions will serve as the constructors for our algebraic expressions. Algebraic expressions are usually formed by operation symbols and then variables as atoms. Instead of variables, our algebraic expressions use terms from our calculus for atoms. The η constructor can thus take an ordinary term from our calculus and make it an atomic algebraic expression. The operation symbols op are then the operations of the algebra.

The other three expression types correspond to functions over algebraic expressions.

-The most useful is the handler . 4 It is an iterator for the type of algebraic expressions. The terms M 1 ,. . . ,M n and M η in op 1 : M 1 , . . . , op n : M n , η: M η are the clauses for the constructors op 1 ,. . . ,op n and η, respectively. We will use handlers to define interpretations of operation symbols in algebraic expressions. -The cherry -• operator allows us to extract terms out of algebraic expressions. If an algebraic expression is of the form η M , applying -• to it will yield M . -The exchange operator C permits a kind of commutation between the λbinder and the operation symbols. We will see its use later.

Types

We now give a syntax for the types of our calculus along with a typing relation.

In the grammar below, α, β, γ. . . range over types, ν ranges over atomic types from some set T and E, E . . . range over effect signatures (introduced below). The types of our language consist of:

α, β, γ ::= α → β [function] | ν [atom] | F E (α) [computation]
The only novelty here is the F E (α) computation type. This is the type of algebraic expressions whose atoms are terms of type α and whose operation symbols come from the effect signature E. We call them computation types and we call terms of these types computations because our algebraic expressions will always represent some kind of program with effects.

Effect signatures are similar to typing contexts. They are partial mappings from the set of operation symbols E to pairs of types. We will write the elements of effect signatures the following wayop : α

β ∈ E means that E maps op to the pair of types α and β. 5 When dealing with effect signatures, we will often make use of the disjoint union operator . The term E 1 E 2 serves as a constraint demanding that the domains of E 1 and E 2 be disjoint and at the same time it denotes the effect signature that is the union of E 1 and E 2 .

The typing rules are presented in Figure 1.

The typing rules mirror the syntax of terms. Again, the first four rules come from STLC. The [η] and [ -• ] rules are self-explanatory and so we will focus on the [op], [ ] and [C] rules.

[op] To use an operation op : α β, we provide the input parameter M p : α and a continuation λx. M c : β → F E (γ), which expects the output of type β. The resulting term has the same type as the body of the continuation, F E (γ).

Before, we have spoken of terms of type F E (γ) as of algebraic expressions generated by the terms of type γ and the operators in the effect signature E. However, having seen the typing rule for operation terms, it might not be obvious how such a term represents an algebraic expression. Traditionally, algebraic signatures map operation symbols to arities, which are natural numbers. Our effect signatures map each operation symbol to a pair of types α β.

-We can explain α by analogy to the single-sorted algebra of vector spaces. In a single-sorted vector space algebra, scalar multiplication is viewed as a unary operation parameterized by some scalar. So technically, there is a different unary operation for each scalar. All of our operations are similarly parameterized and α is the type of that parameter.

Γ, x : α M : β

[abs] Γ λx. M : α → β Γ M : α → β Γ N : α [app] Γ M N : β x : α ∈ Γ [var] Γ x : α c : α ∈ Σ [const] Γ c : α Γ M : α [η] Γ η M : FE(α) Γ Mp : α Γ, x : β Mc : FE(γ) op : α β ∈ E [op] Γ op Mp (λx. Mc) : FE(γ) Γ M : F ∅ (α) [ -• ] Γ -• M : α E = {op i : αi βi}i∈I E f E = E E f [Γ Mi : αi → (βi → F E (δ)) → F E (δ)]i∈I Γ Mη : γ → F E (δ) Γ N : FE(γ) [ ] Γ (op i : Mi)i∈I , η: Mη N : F E (δ) Γ M : α → FE(β) [C] Γ C M : FE(α → β)
Fig. 1: The typing rules for our calculus.

-The type β expresses the arity of the operator. When we say that an operator has arity β, where β is a type, we mean that it takes one operand for every value of β [START_REF] Pretnar | Logic and handling of algebraic effects[END_REF]. We can also think of the operator as taking one operand containing x : β as a free variable.

We can look at the algebraic expression op M p (λx. M c ) as a description of a program that:

interacts with its context by some operator called op to which it provides the input M p and from which it expects to receive an output of type β which it will then bind as the variable x and continue as the program described by M c .

[ ] The banana brackets describe iterators/catamorphisms. 6 In the typing rule, E is the input's signature, E is the output's signature, γ is the input's atom type and δ is the output's atom type. E is decomposed into the operations that our iterator will actually interpret, the other operations form a residual signature E f . The output signature will then still contain the uninterpreted operations E f combined with any operations E that our interpretation might introduce.

[C] We said before that the C function will let us commute λ and operations.

Here we see that, on the type level, this corresponds to commuting the F E ( ) and the α → type constructors.

Reduction Rules

We will now finally give a semantics to our calculus. The semantics will be given in the form of a reduction relation on terms. Even though the point of the calculus is to talk about effects, the reduction semantics will not be based on any fixed evaluation order; any subterm that is a redex can be reduced in any context. The reduction rules are given in Fig. 2.

(λx. M ) N → rule β M [x := N ] λx. M x → rule η M where x / ∈ FV(M ) (op i : Mi)i∈I , η: Mη (η N ) → rule η Mη N (op i : Mi)i∈I
-• (η M ) → rule -• M C (λx. η M ) → rule Cη η (λx. M ) C (λx. op Mp (λy. Mc)) → rule Cop op Mp (λy. C (λx. Mc))
where x / ∈ FV(Mp)

Fig. 2: The reduction rules of our calculus.

We have the β and η rules, which, by no coincidence, are the same rules as the ones found in STLC. The rest are function definitions for , -• and C. By looking at the definition of , we see that it is an iterator. It replaces every occurrence of the constructors op j and η with M j and M η , respectively.

The C function recursively swaps C (λx. ) with op M p (λy. ) using the C op rule. When C finally meets the η constructor, it swaps (λx. ) with η and terminates. Note that the constraint x / ∈ FV(M p ) in rule C op cannot be dismissed by renaming of bound variables. If the parameter M p contains a free occurrence of x, the evaluation of C will get stuck. C is thus a partial function: it is only applicable when none of the operations being commuted with the λ-binder actually depend on the bound variable.

Common Combinators

When demonstrating the calculus in the next section, the following combinators will be helpful. First, we define a sequencing operator. The operator =, called bind, replaces all the α-typed atoms of a F E (α)-typed expression with F E (β)typed expressions. More intuitively, M = N is the program that first runs M to get its result x and then continues as the program N x.

= : F E (α) → (α → F E (β)) → F E (β) M = N = η: N M
The type constructor F E along with the operators η and = form a free monad. Using this monadic structure, we can define the following combinators (variations on application) which we will make heavy use of in Section 3.

• : F E (α → β) → α → F E (β) F • x = F = (λf. η (f x)) • : (α → β) → F E (α) → F E (β) f • X = X = (λx. η (f x)) • : F E (α → β) → F E (α) → F E (β) F • X = F = (λf. X = (λx. η (f x)))
All of these operators associate to the left, so

f • X • Y should be read as (f • X) • Y . Let • : o → o → o
be a binary operator on propositions. We define the following syntax for the same operator lifted to computations of propositions.

• : F E (o) → F E (o) → F E (o) M • N = (λmn. m • n) • M • N
3 Linguistic Phenomena as Effects

Deixis

We will now try to use this calculus to do some semantics. Here is our tectogrammar in an abstract categorial grammar presentation [START_REF] De Groote | Towards abstract categorial grammars[END_REF]. In the semantics for me , we use the speaker operation to retrieve the current speaker and make it available as the value of the variable x. The star ( ) passed to speaker is a dummy value of the unit type 1.

This, and all the semantics we will see in this paper, satisfies a homomorphism condition that whenever M : τ , then M : τ . In our case, N P = F E (ι) and S = F E (o), where ι and o are the types of individuals and propositions, respectively. Of E, we assume that speaker : 1 ι ∈ E, since that is the type of speaker used in our semantics. 7 With this fragment, we can give meanings to trivial sentences like:

(1) John loves Mary.

(2) Mary loves me.

whose meanings we can calculate as:

loves Mary John η (love j m) (1)
loves me Mary speaker (λx. η (love m x))

The meaning of ( 1) is a proposition of type o wrapped in η, i.e. something that we can interpret in a model. As for the meaning of (2), the speaker operator has propagated from the me lexical entry up to the meaning of the whole sentence. We now have an algebraic expression having as operands the propositions love m x for all possible x : ι. In order to get a single proposition which is to be seen as the truth-conditional meaning of the sentence and which can be evaluated in a model, we will need to fix the speaker. We will do so by defining an interpreting handler. withSpeaker : ι → F {speaker:1 ι} E (α) → F E (α) withSpeaker = λsM. speaker: (λxk. k s) M Note that we omitted the η clause in the banana brackets above. In such cases, we say there is a default clause η: (λx. η x).

withSpeaker s loves me Mary η (love m s)

So far, we could have done the same by introducing a constant named me to stand in for the speaker. However, since handlers are part of our object language, we can include them in lexical entries. With this, we can handle phenomena such as direct (quoted) speech, that rebinds the current speaker in a certain scope. Here we elaborated the entry for indirect speech so it is easier to compare with the one for direct speech, which just adds a use of the withSpeaker operator.

(3) John said Mary loves me.

The meaning of sentence (3) depends on the speaker (as testified by the use of the speaker operator) whereas in (4), this dependence has been eliminated due to the use of direct speech.

Quantification

Now we turn our attention to quantificational noun phrases. every, a : N -• N P man, woman :

N every := λN. scope (λc. ∀ • (C (λx. (N • x) → (c x)))) (λx. η x) a := λN. scope (λc. ∃ • (C (λx. (N • x) ∧ (c x)))) (λx. η x) man := η man woman := η woman
The entries for every and a might seem intimidating. However, if we ignore the • , the C, the • and the overline on the logical operator, we get the familiar generalized quantifiers. These decorations are the plumbing that takes care of the proper sequencing of effects.

Note that we make use of the C operator here. In the denotation of a , the term (λx. (N • x) ∧ (c x)) describes the property to which we want to apply the quantifier ∃. However, this term is of type ι → F E (o). In order to apply ∃, we need something of type ι → o. Intuitively, the effects of E correspond to the process of interpretation, the process of arriving at some logical form of the sentence. They should thus be independent of the particular individual that we use as a witness for x when we try to model-check the resulting logical form. This independence allows us use the C operator without fear of getting stuck. Once we arrive at the type F E (ι → o), it is a simple case of using ∃ • to apply the quantifier within the computation type. 89 While the terms that use the scope operator might be complex, the handler that interprets them is as simple as can be. SI = λM. scope: (λck. c k) M Same as with withSpeaker, SI will also be used in lexical items. By interpreting the scope operation in a particular place, we effectively determine the scope of the quantifier. Hence the name of SI, short for Scope Island. If we want to model clause boundaries as scope islands, we can do so by inserting SI in the lexical entries of clause constructors (in our case, the verbs).

loves := λOS. SI ( loves O S)

said is := λCS. SI ( said is C S) said ds := λCS. SI ( said ds C S)
Whenever we use the semantic brackets on the right-hand side of these revised definitions, they stand for the denotations we have assigned previously.

(5) Every man loves a woman.

(6) John said every woman loves me.

(7) John said, "Every woman loves me". loves (a woman) (every man) η (∀x.man x → (∃y.woman y ∧ love x y))

withSpeaker s said is (loves me (every woman)) John η (say j (∀x.woman x → love x s))

said ds (loves me (every woman)) John η (say j (∀x.woman x → love x j))

The calculus offers us flexibility when modelling the semantics. We might choose to relax the constraint that clauses are scope islands by keeping the old 8 Other solutions to this problem include separating the language of logical forms and the metalanguage used in the semantic lexical entries to manipulate logical forms as objects [START_REF] Kiselyov | Applicative abstract categorial grammars[END_REF]. 9 Our C has been inspired by an operator of the same name proposed in [START_REF] De Groote | On Logical Relations and Conservativity[END_REF]: de Groote introduces a structure that specializes applicative functors in a similar direction as monads by introducing the C operator and equipping it with certain laws; our C operator makes the FE type constructor an instance of this structure.

entries for verbs that do not use the SI handler. We might then want to add the SI handler to the lexical entry of said ds , next to the withSpeaker handler, so that quantifiers cannot escape quoted expressions. We might also allow for inverse scope readings by, e.g., providing entries for transitive verbs that evaluate their arguments right-to-left (though then we would have to watch out for crossover effects if we were to add anaphora).

Conventional Implicature

Our goal is to show the modularity of this approach and so we will continue and plug in one more phenomenon into our growing fragment: conventional implicatures, as analyzed by Potts [START_REF] Potts | The logic of conventional implicatures[END_REF]. Specifically, we will focus on nominal appositives.

appos : N P -• N P -• N P best-friend : N P -• N P appos := λXY. X = (λx. SI (η x = Y ) = (λi. implicate i (λz. η x))) best-friend := λX. best-friend • X
In the denotation of the nominal appositive construction, appos, we first evaluate the head noun phrase X : N P to find its referent x : ι. We then want to implicate that x is equal to the referent of Y . The term η x = Y (note the line over =) is the term that computes that referent and gives us the proposition we want. We also want to state that no quantifier from within the appositive Y should escape into the matrix clause and so we wrap this computation in the SI handler to establish a scope island. Finally, we pass this proposition as an argument to implicate and we return x as the referent of the noun phrase.

The point of the implicate operation is to smuggle non-at-issue content outside the scope of logical operators. The contribution of an appositive should survive, e.g., logical negation. 10 The place where we will accommodate the implicated truth-conditions will be determined by the use of the following handler: accommodate :

F {implicate:o 1} E (o) → F E (o) accommodate = λM. implicate: (λik. η i ∧ k ) M
We want conventional implicatures to project out of the common logical operators. However, when we consider direct quotes, we would not like to attribute the implicature made by the quotee to the quoter. We can implement this by inserting the accommodate handler into the lexical entry for direct speech. In [START_REF] Giorgolo | Natural language semantics with enriched meanings[END_REF], the conventional implicature that John is the speaker's best friend projects from the scope of the quantifier. On the other hand, in [START_REF] Hobbs | Making computational sense of montague's intensional logic[END_REF], the implicature does not project from the quoted clause and so it is not misattributed. withSpeaker s (accommodate loves (every woman) (appos John (best-friend me)) ) η ((j = best-friend s) ∧ (∀x.woman x → love j x))

accommodate loves John (appos Mary (best-friend everyone)) η ((∀x.m = best-friend x) ∧ (love m j))

said ds (loves me (appos (best-friend me) Mary)) (a man) η (∃x.man x ∧ say x ((best-friend x = m) ∧ (love (best-friend x) x))) (10)

Summary

Let us look back at the modularity of our approach and count how often during the incremental development of our fragment we either had to modify existing denotations or explicitly mention previous effects in new denotations. When adding quantification:

in the old denotations of verbs, we added the new SI handler so that clauses form scope islands

When adding appositives and their conventional implicatures:

in the old denotations said ds , we added the new accommodate handler to state that conventional implicatures should not project out of quoted speech in the new denotation appos , we used the old SI handler to state that appositives should form scope islands Otherwise, none of the denotations prescribed in our semantic lexicon had to be changed. We did not have to type-raise non-quantificational NP constructors like John , me or best-friend . With the exception of direct speech, we did not have to modify any existing denotations to enable us to collect conventional implicatures from different subconstituents.

Furthermore, all of the modifications we have performed to existing denotations are additions of handlers for new effects. This gives us a strong guarantee that all of the old results are conserved, since applying a handler to a computation which does not use the operations being handled changes nothing.

The goal of our calculus is to enable the creation of semantic lexicons with a high degree of separation of concerns. In this section, we have seen how it can be done for one particular fragment.

The calculus defined in Sect. 2 and to which we will refer as λ , has some satisfying properties.

First of all, the reduction rules preserve types of terms (subject reduction). The reduction relation itself is confluent and, for well-typed terms, it is also terminating. This means that typed λ is strongly normalizing.

The proof of subject reduction is a mechanical proof by induction. For confluence and termination, we employ very similar strategies: we make use of general results and show how they apply to our calculus. Due to space limitations, we will pursue in detail only the proof of confluence.

Our reduction relation is given as a set of rules which map redexes matching some pattern into contracta built up out of the redexes' free variables. However, our language also features binding, and so some of the rules are conditioned on whether or not certain variables occur freely in parts of the redex. Fortunately, such rewriting systems have been thoroughly studied. Klop's Combinatory Reduction Systems (CRSs) [START_REF] Klop | Combinatory reduction systems: introduction and survey[END_REF] is one class of such rewriting systems.

We will make use of the result that orthogonal CRSs are confluent [START_REF] Klop | Combinatory reduction systems: introduction and survey[END_REF]. A CRS is orthogonal if it is left-linear and non-ambiguous. We will need to adapt our formulation of the reduction rules so that they form a CRS and we will need to check whether we satisfy left-linearity and non-ambiguity (we will see what these two properties mean when we get to them).

We refer the reader to [START_REF] Klop | Combinatory reduction systems: introduction and survey[END_REF] for the definition of CRSs. The key point is that in a CRS, the free variables which appear in the left-hand side of a rewrite rule, called metavariables, are explicitly annotated with the set of all free variables that are allowed to occur within a term which would instantiate them. This allows us to encode all of our x / ∈ F V (M ) constraints. One detail which must be taken care of is the set notation (op i : M i ) i∈I and the indices I used in the rules. We do away with this notation by adding a separate rule for every possible instantiation of the schema. This means that for each sequence of distinct operation symbols op 1 ,. . . ,op n , we end up with:

a special rewriting rule op

1 : M 1 , . . . , op n : M n , η: M η (η N ) → M η N -for every 1 ≤ i ≤ n, a special rewriting rule op 1 : M 1 , . . . , op n : M n , η: M η (op i N p (λx. N c (x))) → M i N p (λx. op 1 : M 1 , . . . , op n : M n , η: M η N c (x)) -for every op ∈ E \ {op i |1 ≤ i ≤ n}, a special rewriting rule op 1 : M 1 , . . . , op n : M n , η: M η (op N p (λx. N c (x))) → op N p (λx. op 1 : M 1 , . . . , op n : M n , η: M η N c (x))
So now we have a CRS which defines the same reduction relation as the rules we have shown in 2.3. Next, we verify the two conditions. Left-linearity states that no left-hand side of any rule contains multiple occurrences of the same metavariable. By examining our rules, we find that this is indeed the case. 11 Non-ambiguity demands that there is no non-trivial overlap between any of the left-hand sides. 12 In our CRS, we have overlaps between the β and the η rules. We split our CRS into one with just the η rule (→ η ) and one with all the other rules (→ λ ). Now, there is no overlap in either of these CRSs, so they are both orthogonal and therefore confluent.

We then use the Lemma of Hindley-Rosen [17, p. 7] to show that the union of → λ and → η is confluent when → λ and → η are both confluent and commute together. For that, all that is left to prove is that → λ and → η commute. Thanks to another result due to Hindley [17, p. 8], it is enough to prove that for all a, b and c such that a → λ b and a → η c, we have a d such that b η d and c → = λ d. The proof of this is a straightforward induction on the structure of a.

Conclusion

In our contribution, we have introduced a new calculus motivated by modelling detailed semantics and inspired by current work in programming language theory. Our calculus is an extension of the simply-typed lambda calculus which is the de facto lingua franca of semanticists. Its purpose is to facilitate the communication of semantic ideas without depending on complex programming languages [START_REF] Maršík | Algebraic Effects and Handlers in Natural Language Interpretation[END_REF][START_REF] Kiselyov | Lambda: the ultimate syntax-semantics interface[END_REF] and to do so with a well-defined formal semantics.

We have demonstrated the features of our calculus on several examples exhibiting phenomena such as deixis, quantification and conventional implicature. While our calculus still requires us to do some uninteresting plumbing to be able to correctly connect all the denotations together, we have seen that the resulting denotations are very generic. We were able to add new phenomena without having to change much of what we have done before and the changes we have made arguably corresponded to places where the different phenomena interact.

Finally, we have also shown that the calculus shares some of the useful properties of the simply-typed lambda calculus, namely strong normalization.

In future work, it would be useful to automate some of the routine plumbing that we have to do in our terms. It will also be important to test the methodology on larger and more diverse fragments (besides this fragment, we have also created one combining anaphora, quantification and presupposition [START_REF] Maršík | Algebraic Effects and Handlers in Natural Language Interpretation[END_REF]). Last but not least, it would be interesting to delve deeper into the foundational differences between the approach used here, the monad transformers used by Charlow [START_REF] Charlow | On the semantics of exceptional scope[END_REF] and the applicative functors used by Kiselyov [START_REF] Kiselyov | Applicative abstract categorial grammars[END_REF].

  , η: Mη (op j Np (λx. Nc)) → rule op Mj Np (λx. (op i : Mi)i∈I , η: Mη Nc) where j ∈ I and x / ∈ FV((Mi)i∈I , Mη) (op i : Mi)i∈I , η: Mη (op j Np (λx. Nc)) → rule op op j Np (λx. (op i : Mi)i∈I , η: Mη Nc) where j / ∈ I and x / ∈ FV((Mi)i∈I , Mη)

  John, Mary, me : N P loves : N P -• N P -• S And here is our semantics. John := η j Mary := η m me := speaker (λx. η x) loves := λOS. love • S • O

  said is : S -• N P -• S said ds : S -• N P -• S Those are our new syntactic constructors: one for the indirect speech use of said and the other for the direct speech use (their surface realizations would differ typographically or phonologically). Let us give them some semantics. said is = λCS. say • S • C = λCS. S = (λs. say s • C) said ds = λCS. S = (λs. say s • (withSpeaker s C))

( 4 )

 4 John said, "Mary loves me". said is (loves me Mary) John speaker (λx. η (say j (love m x))) (3) said ds (loves me Mary) John η (say j (love m j))

  said ds := λCS. SI (S = (λs. say s • (withSpeaker s (accommodate C)))) Consider the following three examples. (8) John, my best friend, loves every woman. (9) Mary, everyone's best friend, loves John.(10) A man said, "My best friend, Mary, loves me".

Side effects are to programming languages what pragmatics are to natural languages: they both study how expressions interact with the worlds of their users. It might then come as no surprise that phenomena such as anaphora, presupposition, deixis and conventional implicature yield a monadic description.

Also known as monad transformers in functional programming.

Pronounced "banana". See[START_REF] Meijer | Functional programming with bananas, lenses, envelopes and barbed wire[END_REF] for the introduction of banana brackets.

The two types α and β are to be seen as the operation's input and output types, respectively.

These are similar to recursors/paramorphisms. See[START_REF] Meijer | Functional programming with bananas, lenses, envelopes and barbed wire[END_REF] for the difference. Catamorphisms are also known as folds and the common higher-order function fold found in functional programming languages is actually the iterator/catamorphism for lists.

The definition of (non-trivial) overlap is the same one as the one used when defining critical pairs. See[START_REF] Klop | Combinatory reduction systems: introduction and survey[END_REF] for the precise definition.


Nicolas Cordier

Hervé Delingette

Matthieu Lê

Nicholas Ayache

Extended Modality Propagation: Image Synthesis of Pathological Cases

Keywords: patch-based, multi-atlas, glioma, generative model, medical image simulation, modality synthesis

This paper describes a novel generative model for the synthesis of multi-modal medical images of pathological cases based on a single label map. Our model builds upon i) a generative model commonly used for label fusion and multi-atlas patch-based segmentation of healthy anatomical structures, ii) the Modality Propagation iterative strategy used for a spatiallycoherent synthesis of subject-specific scans of desired image modalities. The expression Extended Modality Propagation is coined to refer to the extension of Modality Propagation to the synthesis of images of pathological cases. Moreover, image synthesis uncertainty is estimated. An application to Magnetic Resonance Imaging synthesis of glioma-bearing brains is i) validated on the training dataset of a Multimodal Brain Tumor Image Segmentation challenge, ii) compared to the state-of-theart in glioma image synthesis, and iii) illustrated using the output of two different tumor growth models. Such a generative model allows the generation of a large dataset of synthetic cases, which could prove useful for the training, validation, or benchmarking of image processing algorithms.

I. INTRODUCTION

A. Motivation

The availability of public datasets [START_REF] Jack | The Alzheimer's disease neuroimaging initiative (ADNI): MRI methods[END_REF]- [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF] of annotated medical images is a key factor in the development of medical image computing. For instance, the organization of the Multimodal Brain Tumor Image Segmentation (BraTS) benchmark challenge [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF] has lead to an important and fruitful research activity in glioma segmentation. However, the creation of a benchmark dataset is costly for obvious reasons: i) a large number of cases is required in the training and testing datasets to capture the variability of structures or pathologies, ii) images should be annotated by experts, which requires a complex and time-consuming manual work, and can still lead to debatable results (e.g. inter-rater variability in the range 74-85% for glioblastoma [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF], datasets can contain incorrect segmentations [START_REF] Wachinger | Brainprint: A discriminative characterization of brain morphology[END_REF]), iii) the distribution of medical images leads to ethical concerns, since the identification of patients may be possible despite anonymization steps such as skullstripping [START_REF] Wachinger | Brainprint: A discriminative characterization of brain morphology[END_REF], iv) and finally, ensuring the continued quality of a very large dataset of anonymized images is complex, due to the presence of longitudinal data and duplicates [START_REF] Toews | A featurebased approach to big data analysis of medical images[END_REF]. The development of image synthesis could allow i) the augmentation of a dataset by including new realistic synthetic cases, and ii) the Université Côte d'Azur, Inria, France. Asclepios Research Project ; email address: herve.delingette@inria.fr. creation of large, annotated, unbiased, anonymized, and easyto-maintain datasets of synthetic medical images of virtual patients (namely atlases completely generated by a computational biophysical model). The availability of a large dataset is especially important for pathologies such as glioblastoma which exhibit a high variability of shape and appearance. For reference, the 2013 BraTS challenge only consisted of 20 real high-grade cases for training, and 10 real high-grade cases for benchmarking [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF]). The 2014 BraTS challenge included about 250 additional cases obtained from The Cancer Imaging Archive [START_REF] Clark | The cancer imaging archive (tcia): maintaining and operating a public information repository[END_REF], however other problems arose: the ground truth was unavailable and had to be obtained through a consensus of automatic glioma segmentation algorithms, which resulted in a ground truth of lower reliability [START_REF] Havaei | Brain tumor segmentation with deep neural networks[END_REF].

B. Related work

The problem which is tackled in this paper is the synthesis of multi-modal medical images of pathological cases, based on a single label map, as illustrated in Figure 1 for brain tumors imaged with Magnetic Resonance Imaging (MRI). In the following, the focus is on related work regarding MRI synthesis of healthy and tumor-bearing brains.

1) MRI synthesis of healthy brains: Original attempts at generating synthetic MRI of healthy brains relied on a numerical simulation of MR acquisition physics. Given MR scanner parameters such as echo time and relaxation time, a discreteevent simulation model [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF]- [START_REF] Glatard | A Virtual Imaging Platform for multi-modality medical image simulation[END_REF] describes the dynamics of magnetization vectors, at each spatial position, according to Bloch equations [START_REF] Bloch | Nuclear induction[END_REF]. Such numerical simulations are computationally expensive, except for specific cases for which there exists a closed-form solution to Bloch equations, such as for spin-echo or gradient-echo MR sequences [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF], [START_REF] Maitra | Synthetic magnetic resonance imaging revisited[END_REF], [START_REF] Iglesias | Synthetic MRI signal standardization: application to multi-atlas analysis[END_REF].

To describe the imaged object, two strategies are possible. The first strategy consists in describing the geometry of the imaged tissues with tissue-specific templates, and providing the tissue-specific biophysical properties (spin density and relaxation times). This allows the synthesis of MRI acquired on virtual geometrical templates. However, it requires the knowledge of tissue-specific biophysical properties [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF], which are poorly referenced in the literature for glioblastoma compartments. The second strategy does not rely on the definition of tissues: biophysical properties are specified voxel-wise, after an estimation from several MR scans obtained in a short time-frame with a very strict acquisition protocol (quantitative MRI or relaxometry), or by optimization methods [START_REF] Maitra | Synthetic magnetic resonance imaging revisited[END_REF], [START_REF] Iglesias | Synthetic MRI signal standardization: application to multi-atlas analysis[END_REF]. This strategy does not allow to generate synthetic images on virtual patients. Moreover, it cannot be applied to tumor scans when a proton density map is not available.

Recently, an iterative patch-based heuristic [START_REF] Ye | Modality Propagation: coherent synthesis of subject-specific scans with data-driven regularization[END_REF] was used for Modality Propagation, i.e. the synthesis of a realistic subjectspecific scan of one modality given a scan acquired with another modality. Roots of the Modality Propagation framework include image reconstruction works, e.g. PET image reconstruction incorporating MR data as prior information for anatomical boundary [START_REF] Leahy | Incorporation of anatomical MR data for improved functional imaging with PET[END_REF] and example-based approaches using non-local patch similarity constraints, known as brain image hallucination [START_REF] Rousseau | Brain hallucination[END_REF]- [START_REF] Roy | MR contrast synthesis for lesion segmentation[END_REF]. The core of the algorithm consists in modeling the covariation of local intensity patterns across modalities, in a fashion inspired by multi-atlas patchbased segmentation of healthy brain structures [START_REF] Rousseau | A supervised patchbased approach for human brain labeling[END_REF], [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF]. Successful applications of Modality Propagation or related methods include synthesis of a patient-specific attenuation map for hybrid MR-PET [START_REF] Burgos | Attenuation correction synthesis for hybrid PET-MR scanners: application to brain studies[END_REF], more accurate registration [START_REF] Iglesias | Is synthesizing MRI contrast useful for inter-modality analysis?[END_REF], super-resolution [START_REF] Roy | Magnetic resonance image examplebased contrast synthesis[END_REF], and outlier detection [START_REF] Ye | Modality Propagation: coherent synthesis of subject-specific scans with data-driven regularization[END_REF], [START_REF] Roy | Magnetic resonance image examplebased contrast synthesis[END_REF], [START_REF] Cardoso | Template-Based Multimodal Joint Generative Model of Brain Data[END_REF]. However, Modality Propagation shows two limitations for the synthesis of MRI for virtual pathological cases: i) the local search window framework is adopted, which does not take into account the variability of lesion position, and might restrict the method to the image synthesis of healthy cases, ii) the input is a real image, which is assumed not to be available for virtual cases.

2) MRI synthesis of tumor-bearing brains: Previous works regarding MRI synthesis of tumor-bearing brains [START_REF] Rexilius | A framework for the generation of realistic brain tumor phantoms and applications[END_REF]- [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF] build upon tissue-specific templates, as in the framework introduced for healthy brains [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF]. Typically, a tumor seed is artificially positioned in a healthy brain atlas, then a tumor growth model simulates the tumor extension over time and its mass effect, i.e. the displacement of neighboring healthy tissues, which results in templates of healthy tissues and tumor compartments. However, in contrast to [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF], tissue-specific biophysical properties are not considered in favor of direct specifications of tissue-specific average MR intensities [START_REF] Rexilius | A framework for the generation of realistic brain tumor phantoms and applications[END_REF], [START_REF] Prastawa | Synthetic ground truth for validation of brain tumor MRI segmentation[END_REF]. In [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF], textures are learnt for each healthy tissue and tumor compartment, which results in a set of 3D texture images. The only difference between different synthetic MRI lies in tissue-specific templates: tissue-specific average MR intensities and 3D texture images are therefore always identical.

This last approach is applicable to a wide variety of cases, but also bears some limitations: i) the complex distribution of image intensities for tumor compartments is summarized by its expectation, which is oblivious of multi-modal intensity distributions, ii) the inter-patient MR normalization procedure is not specified, which makes it difficult to standardize real MRI so that they look like synthetic MRI, typically for the training of machine learning algorithms [START_REF] Geremia | Brain tumor cell density estimation from multi-modal MR images based on a synthetic tumor growth model[END_REF], iii) simulated images do not show the variability of intensity of realistic MR scans, and the addition of a very high Gaussian noise only limits this effect.

The simulator of synthetic pathological MRI1 described in [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF] has been used in a number of research articles mostly for prototyping and validation, in the context of glioma segmentation [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF], [START_REF] Geremia | Brain tumor cell density estimation from multi-modal MR images based on a synthetic tumor growth model[END_REF], [START_REF] Hamamci | Tumorcut: segmentation of brain tumors on contrast enhanced MR images for radiosurgery applications[END_REF], outlier detection algorithm [START_REF] Galimzianova | Robust estimation of unbalanced mixture models on samples with outliers[END_REF], registration of a healthy brain atlas to a tumor-bearing patient image [START_REF] Bauer | Multiscale modeling for image analysis of brain tumor studies[END_REF], and construction of a brain atlas [START_REF] Liu | Low-rank Atlas Image Analyses in the Presence of Pathologies[END_REF]. Other applications include the training of machine-learning algorithms for glioma segmentation. [START_REF] Geremia | Brain tumor cell density estimation from multi-modal MR images based on a synthetic tumor growth model[END_REF] trained random forest on a large dataset consisting of 740 synthetic cases and showed good segmentation results on a few real cases for testing. [START_REF] Galimzianova | Robust estimation of unbalanced mixture models on samples with outliers[END_REF] performed a thorough validation of an outlier detection algorithm, based on 100 MR scans synthesized with different tumor volumes to test the robustness of the algorithm to the amount of outliers. Image synthesis could also be useful in the context of tumor growth modeling: in [START_REF] Gu | Applying a patient-specific bio-mathematical model of glioma growth to develop virtual [18F]-FMISO-PET images[END_REF], PET images are synthesized for tumor-bearing brains, so as to allow clinically relevant interpretations of tumor growth model outputs. Although the applications of MRI synthesis of pathological cases are numerous, the competitors in the BraTS benchmark challenge tend to completely ignore the provided synthetic cases [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF], which is likely due to their lack of realism.

In this paper, we address the image synthesis of pathological cases by using an iterative multi-atlas patch-based algorithm, inspired by recent successful algorithms in MRI synthesis of healthy brains. The expression Extended Modality Propagation (EMP) is coined for two reasons: i) to refer to the extension of Modality Propagation for the synthesis of images of pathological cases, ii) in contrast to Modality Propagation, EMP can deal with label maps in addition to image intensity. The synthesis process only requires a single label map, or the output of any tumor growth model in terms of the usual tumor compartments. Realistic synthetic images are obtained in the specific application of our algorithm to MRI synthesis of tumor-bearing brains. Our contributions include a novel generative model, a heuristic iterative algorithm to solve for the posterior distribution of multi-channel MR intensities, the estimation of image synthesis uncertainty, the analysis of Modality Propagation iterative feature augmentation, and the public availability of a large dataset of annotated and realistic MRI exhibiting gliomas. In the remainder, the generative model is described and solved for the maximum a posteriori and for uncertainty estimation (Section II). Then synthesized MRI are compared to real MRI from the training dataset of BraTS benchmark challenge, and to the state-of-the-art in glioma image synthesis, and illustrated using the output of two different tumor growth models (Section III). Finally, results and perspectives are discussed (Sections IV and V).

II. EXTENDED MODALITY PROPAGATION

In the following, a label map L describing the anatomy of a tumor-bearing brain is assumed to be known. The goal of Fig. 2: Pipeline for Extended Modality Propagation. Green disks mark known tumor locations. A fundamental idea is that central voxels of similar label patches should have similar multi-channel image intensity. During the first iteration, multichannel intensity J is assigned for every test label patch based on comparisons with a database of training label patches. During subsequent iterations, multi-channel intensity J is updated based on comparisons between augmented label patches, which are the concatenation of label patches and multi-channel image patches synthesized during the previous iteration (feedback).

Extended Modality Propagation (EMP) is to synthesize medical images corresponding to the same anatomy. To achieve this goal, EMP relies on a set of N training cases {L n , I n } n where L n is a label map and I n denotes multi-channel MRI (T1, T1C, T2, FLAIR). In practice, the set of labels comprises L = 7 classes: cerebrospinal fluid, grey matter, white matter, necrotic core, edema, non-enhancing tumor core, and enhancing core.

The EMP probabilistic generative model is first introduced to describe the synthesis of images. Second, an approximation of the marginalization step is introduced to estimate the maximum a posteriori and the image synthesis model uncertainty. Finally, more consistent patch matching is obtained by augmenting label patches: i) multi-scale patches allow to consider larger patches in a computationally efficient manner, ii) a patch match invariant to cube isometries allows to consider larger patches without being impacted by a smaller sample size, and iii) as for Modality Propagation, an iterative patch match procedure results in smoother and more consistent synthetic images. The EMP pipeline is detailed in Figure 2.

A. Definition of a probabilistic generative model 1) Graphical model: The proposed generative model for Extended Modality Propagation builds upon the Bayesian modeling of local weighted voting for multi-atlas segmentation [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF], which was then adapted for patch-based segmentation [START_REF] Bai | A probabilistic patch-based label fusion model for multi-atlas segmentation with registration refinement: application to cardiac MR images[END_REF], [START_REF] Iglesias | Multi-Atlas Segmentation of Biomedical Images: A Survey[END_REF]. Conventionally, a patch S (I, x) is a cube of edge length 2r + 1, centered at spatial position x ∈ R 3 , consisting of image intensity values taken from an image I.

Our model is based on the 5 following assumptions. i) Every patch S (L, x) originates from a single patch S (L n , y) from the training database. This means that the target segmentation L is the result of sampling patches from atlases (L n , I n ) at different positions y in the reference space Ω. The membership index M encodes, for every position x ∈ Ω, both the atlas index n and the position y to sample from: M(x) = (n, y). ii) There is a uniform prior on M: no training case n or patch location y is preferred a priori. This assumption is in sharp contrast with a common assumption of local search window in the patch-based framework. iii) The only constraint on M is that the source patch for a target patch S (L, x) has the same central label L n (y) = L(x). iv) If M were known, the target image J would be given by intensity, perturbed with white noise, I n (y) at central voxels of sampled patches. v) Given M(x), probabilities of observing intensity J(x) and label L(x) are conditionally independent. This assumption introduces a symmetry between segmentation and synthesis problems. Moreover, this assumption is reasonable: if M(x) were known, then the information regarding label would be completely defined by L n only, and the information regarding intensity would be completely defined by I n only.

The corresponding graphical model is shown in Figure 3. This model mirrors the graphical model described in [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF] for multi-atlas patch-based segmentation of pathological cases.

2) Class-conditional image intensity likelihood: Let j ∈ R d denote multi-channel voxel-wise image intensity, with d being the number of MR channels. In our case, d = 4. Class-conditional image intensity likelihood p n J(x) = j | I n , y, Σ Ln(y) , L n (y), L(x) is the probability that the target intensity J(x) = j is sampled from a multivariate Gaussian probability density function N , with

M L J L n {I n } σ 2 Σ l

|Ω|

Fig. 3: Graphical model for Extended Modality Propagation.

A membership variable M : Ω → {1, . . . , N } × Ω is sampled at every position x in the MNI space Ω to encode the training case n to sample from, and the spatial offset y: M(x) = (n, y). The observed Log-Odds map L(x) results from the sampling of a patch S( L, x) from a Gaussian distribution with S( L n , y) as mean and isotropic scaling σ 2 (x) as variance. The intensity J at position x is sampled from a Gaussian distribution with the central value I n (y) as mean and Σ L(x) as covariance matrix.

I n (y) as mean and Σ Ln(y) = Σ L(x) as covariance matrix:

p n J(x) = j | I n , y, Σ Ln(y) , L n (y), L(x) ∝ δ Ln(y),L(x) × N j | I n (y), Σ Ln(y)
where δ •,• is Kronecker delta, which translates assumption iii.

In this case, L n (y) = L(x), and the notation for intensity likelihood is simplified to p n J(x) = j | I n , y, Σ L(x) .

3) Log-Odds label maps: The likelihood label model relies on continuous Log-Odds maps L and L n [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF], i.e. probability values computed from categorical label maps L and L n . Indeed, during the fist iteration of image synthesis, patch match only relies on a label map as input. It is crucial to choose the best representation for this single input. A categorical label map could be used, but this representation does not distinguish patches inside a tissue from those at its border. If categorical label maps were directly used, the matching of similar label configurations could be poor for label configurations not present in the training dataset. A better strategy, as recommended in [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF], consists in adopting Log-Odds, so as to take into account the distance to borders between tissues and to be able to distinguish different patches among patches completely enclosed in a tissue. Log-Odds distances are defined as proportional to exp ρ D l (x) , where ρ > 0 is the slope constant (in our experiments, ρ = 1 as in [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF]), and D l (L, x) is the signed distance transform of label l for label map L at spatial position x, which is positive inside structures labelled as l.

4) Augmented Log-Odds label maps: During the first iteration t = 1 of the algorithm, the notation L corresponds to a multi-label Log-Odds map: a vector in R L consisting of Log-Odds label probability values is assigned for each spatial position x. Therefore, a patch of edge length r would result in a vector in R (2r+1) 3 ×L . However, for subsequent iterations t ≥ 2, the same notation L corresponds to augmented Log-Odds maps: for each spatial position x, the vector in R L of Log-Odds label probability values is concatenated with the vector in R d of synthetic image intensity of the previous iteration. In this case, a patch of edge length r would result in a vector in R (2r+1) 3 ×(L+d) .

5) Label patch likelihood: Label patch likelihood

p n S( L, x) | L n , y, σ 2 (x)
comprises the similarity between label patches. It is defined as the probability that the target Log-Odds patch S( L, x) is sampled from a multivariate Gaussian probability density function, with S( L n , y) as mean and isotropic scaling σ 2 (x) as variance:

p n S( L, x) | L n , y, σ 2 (x) ∝ N S( L, x) | S( L n , y), σ 2 (x) × 1
where 1 is the identity matrix.

6) Symmetry:

The symmetry between the segmentation model described in [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF] and our synthesis model is emphasized by the definitions of label and intensity likelihoods, as shown in Table I. Indeed, regarding likelihood for observed variables, image intensity likelihood for segmentation [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF] and label patch likelihood for synthesis are Gaussian densities on patches, respectively image patches and Log-Odds patches, with an isotropic scaling parameter for variance, respectively σ 2 n (x) and σ 2 (x). Moreover, regarding likelihood for hidden variables, label likelihood for segmentation [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF] and classconditional image intensity likelihood for synthesis are defined voxel-wise instead of patch-wise, and both involve a Kronecker delta.

B. Inference

Since M is not observed, the inference has to be performed by marginalizing over M [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF]. The marginalization consists in a weighted vote involving patches at every position y in every atlas n, with weights proportional to both label patch and intensity likelihoods.

Variance parameters are first estimated ; second, the problem of estimating J is stated by marginalizing over M ; third, a strategy to determine the arg-max of the posterior distribution is proposed ; finally, a method to estimate the uncertainty of the synthesis procedure is mentioned.

1) Estimation of the variance parameters: There are two kinds of variance parameters: i) covariance matrices {Σ l } l∈{1, ..., L} for class-conditional image intensity likelihoods, which are pre-computed, and ii) isotropic scaling σ 2 (x) for label patch likelihood, which is computed online.

For each label l ∈ {1, . . . , L}, a covariance matrix Σ l results from the fit, using Expectation-Maximization [START_REF] Liu | ML estimation of the t distribution using EM and its extensions, ECM and ECME[END_REF], of a multivariate Student distribution to average intensity of pure patches, i.e. for which all voxels belong to the same class. Robustness to outliers is achieved with Student distributions to model the heavy-tailed distributions of multi-channel image intensity. Supplementary material, available in the supplementary files /multimedia tab, presents 2D projections of obtained Student distributions.

As with patch-based segmentation, the variance σ 2 (x) for Log-Odds patch is spatially varying. In contrast to [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF], σ 2 (x)

Likelihoods

Segmentation model Synthesis model

Intensity likelihood N S(J, x) | S(I n , y), σ 2 n (x) × 1 δ Ln(y),L(x) × N J (x) | I n (y), Σ Ln(y) Label likelihood δ Ln(y), L(x) N S( L, x) | S( L n , y), σ 2 (x) × 1
is a heuristic inspired by [START_REF] Coupé | Patch-based segmentation using expert priors: Application to hippocampus and ventricle segmentation[END_REF], [START_REF] Wang | Geodesic patch-based segmentation[END_REF] and involving the minimal patch distance, defined such that:

√ 2σ (x) = min n,y d S( L, x), S( L n , y)
The same formula is applied to augmented Log-Odds patches during subsequent iterations of the proposed algorithm. As proximity between samples is not well characterized by their distance in high dimensional spaces, the use of the minimal distance instead of another measure (such as the average distance) is preferred to limit the number of patches with high probability.

2) Approximate marginalization: Assuming the membership index M(x) is independent and identically uniformly distributed p (M) = (N × |Ω|)

-|Ω| , the marginalization over M [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF] is given voxel-wise by:

J(x) = argmax j N n=1 y∈Ω p n S( L, x) | L n , y, σ 2 (x) ×p n J(x) = j | I n , y, Σ Ln(y) , L n (y), L(x)
The marginalization over M would require a high number of comparisons between patches, which is infeasible in practice. However, the marginalization step could be approximated, in a similar fashion as [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF]. Let y n (S(L, x)) be the position in reference space of the closest Log-Odds patch, in canonical L 2 distance, of central label l = L(x) found in the atlas (L n , I n ). For image synthesis of pathological cases, we propose to marginalize over only N patches, where N is the total number of atlases: given a target Log-Odds patch of central label L(x), the approximate marginalization only includes patches found at the spatial position y n (S(L, x)) in atlas n, which results in:

J(x) = argmax j N n=1 p n S( L, x) | L n , y n (S(L, x)) , σ 2 (x) ×p n J(x) = j | I n , y n (S(L, x)) , Σ L(x)
3) Arg-max of the posterior distribution: With the approximate marginalization, our original belief was that J(x) could be efficiently computed thanks to Banach Fixed Point Theorem. Indeed, a necessary condition for j * to be an optimal solution is given by gradient cancelling. Let F : R 4 → R be:

F (j) = N n=1 w n × p n J(x) = j | I n , y n (S(L, x)) , Σ L(x)
where:

w n = p n S( L, x) | L n , y n (S(L, x)) , σ 2 (x) > 0
Gradient cancelling gives the necessary condition that j * is a fixed point of a function T : R 4 → R 4 :

j * = T (j * ) = N n=1 λ n (j * ) × µ n
where:

µ n = I n (y n (S(L, x))) λ n (j) = w n × p n J(x) = j | I n , y n (S(L, x)) , Σ L(x) N n=1 w n × p n J(x) = j | I n , y n (S(L, x)) , Σ L(x) ∀n, λ n (j) ≥ 0 N n=1 λ n (j) = 1
For voxels x for which T is a contraction mapping, Banach Fixed Point Theorem gives that j * is unique, and that for any initialization j 0 ∈ R 4 , any sequence (j k ) k∈N such that j k+1 = T (j k ) converges to the fixed point j * with a geometric rate. An implementation of the EMP algorithm is given in Algorithm 1. In fact, we adopt a slight variant of this approach where λ n (j) is computed with a more robust Student likelihood (Section II-B1) instead of a Gaussian likelihood p n J(x) = j | I n , y n (S(L, x)) , Σ L(x) , to cope with the rare case where a few µ n are far away from j * . We could not prove that the Lipschitz constant of T is strictly less than 1 for every voxel in the image, hence the assumption of Banach Fixed Point Theorem may not be fulfilled. Nonetheless, in our experiments for all test cases, the algorithm does converge in about K = 60 iterations for 99.9% of the voxels of the brain mask. For the remaining 0.01% of the voxels, the intensity value obtained at the last iteration is used.

4) Estimation of the image synthesis uncertainty: With the Fixed Point procedure, the estimation of the arg-max of the posterior image intensity distribution involves the iterative estimation of j * as a weighted average of intensities µ n found in the atlases. Therefore, the uncertainty of the image synthesis process can be similarly estimated as the weighted sample variance σ 2 , thanks to the optimal weights λ n (j * ). Estimation of image synthesis uncertainty is illustrated in Figure 4. Uncertainty estimation could allow to pin-point regions with unlikely image intensity based on the proposed model.

C. Consistent patch match

Patch match consistency can be improved by relying on augmented Log-Odds patches and by adopting an iterative feature augmentation strategy.

1) Multi-scale patches: As shown in [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF], multi-scale patches are computationally efficient representations for large patches, which allows more effective patch match. The central part of the patch is described voxel-wise, while the peripheral part is described by average intensity values over neighboring patches, by analogy with the foveal vision. Similar ideas regarding foveation and non-uniform sampling are detailed in [START_REF] Ciresan | Deep neural networks segment neuronal membranes in electron microscopy images[END_REF]. In the following, parameters are set as in [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF]: a 3x3x3 patch is used for the central part, and a 26-neighborhood for the peripheral part.

2) Patch match with invariance properties: Every atlas is affinely registered to the MNI space [START_REF] Mazziotta | A probabilistic atlas and reference system for the human brain: International Consortium for Brain Mapping (ICBM)[END_REF] so that patch extraction is invariant to the pose of the subject. Moreover, the canonical L 2 distance used for patch matching is sensitive to rotation or symmetry of the patches. As shown in [START_REF] Cordier | A patch-based approach for the segmentation of pathologies: Application to glioma labelling[END_REF],

the combination of multi-scale patches and invariance to cube isometries yields better segmentation results, and thus could be beneficial for patch-based image synthesis as well. The 48 cube isometries consist of certain rotations under which the cube is invariant, plus their composition with central symmetry. For healthy brain structures, rotated patches can be unrealistic, only sagittal plane symmetry is considered. For pathological patches, the similarity measure should be insensitive to rotation and symmetry: the number of atlases is small, and using multi-scale patches further decreases the sample size. Therefore, the addition of plausible configurations of patches, generated by applying sagittal plane symmetry or cube isometries, to observed patches is relevant, and ultimately leads to a patch match invariant to the group of cube isometries. Indeed, if we were to apply cube isometries to the target patch, the minimal patch distance between the transformed target patch and the augmented learning dataset would be equal to the minimal patch distance between the target patch and the augmented learning dataset.

3) Iterative feature augmentation: An iterative augmentation procedure guarantees a spatially more coherent patch match, which results in an overall more coherent image synthesis. This is illustrated in Figure 5. Augmented Log-Odds patches are the concatenation of i) Log-Odds patches, weighted with 1 -α t , and ii) multi-channel image patches synthesized during the previous iteration t -1, weighted with α t . There are different possible definitions for weight α t . As with Modality Propagation [START_REF] Ye | Modality Propagation: coherent synthesis of subject-specific scans with data-driven regularization[END_REF], the feedback weight α t increases with the number of iterations, starting with α 1 = 0 for iteration t = 1. Unless mentioned, a maximal number of iterations is set to t max = 3, and α t = (t -1) /t max , so that the weight of Log-Odds patches is never zero and decreases linearly with the number of iterations.

III. RESULTS

Our goal is to synthesize four MR channels commonly acquired for glioma assessment: pre-contrast T 1 -weighted image (T1), contrast-enhanced T 1 -weighted (T1C), T 2 -weighted (T2), and T 2 -FLAIR MR images. Fast approximate nearestneighbour search are used for patch match, more precisely multiple randomized k-d trees for high dimensional data [START_REF] Muja | Scalable Nearest Neighbor Algorithms for High Dimensional Data[END_REF].

A. Pre-processing pipeline

Every image is affinely registered to an MNI atlas and warped to the same reference space. For experiments which require strictly more than 3 iterations, a resampling to 2 mm isotropic resolution is performed to decrease computation time, otherwise we keep the original 1 mm isotropic resolution. A resampling to 2 mm resolution was used by the top-performing glioma segmentation method at the 2012 BraTS challenge, which indicates that the generation of a synthetic dataset of 2 mm resolution images could already be of interest for applications other than medical image synthesis. Atlas images I n are standardized with a global affine intensity transform, applied to each modality independently: non-zero intensity are clipped below 1% and above 99% quantiles, then intensity mean and standard deviation are set to the corresponding Fig. 5: From left to right: input label map, synthetic contrast-enhanced T 1 -weighted MRI after iterations t = 1, 2, 3, and ground truth (real MRI). Image parts near the brain border are better synthesized due to iterative feature augmentation. The presence of dark synthesized regions around the brain, instead of the presence of skull, is due to the fact that atlases are skullstripped. Supplementary material, available in the supplementary files /multimedia tab, presents maps of absolute differences between synthesized images and ground truth. average values over the training dataset, respectively 360 as mean and 120 as standard deviation. Except for the clipping, and potentially the resampling if performed, all the proposed pre-processing operations can be reverted, which eases the adaptation of such a synthetic dataset for the analysis of another real dataset.

B. Validation & comparison to state-of-the-art 1) Dataset and evaluation:

The dataset consists of MRI of 20 high-grade (HG) and 10 low-grade (LG) cases, skullstripped and made publicly available as training dataset in the 2013 edition of the MICCAI BraTS benchmark [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF]. A leave-one-out procedure is always applied to exclude the target image from the training dataset. There is no atlas pre-selection and all training cases, remaining after leave-one-out, are used as atlases for the synthesis of MRI of cases of the same grade. Only the segmentation of the tumor compartments is originally available [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BraTS)[END_REF], so healthy tissues are automatically segmented by applying FSL FAST [START_REF] Zhang | Segmentation of brain MR images through a hidden Markov random field model and the expectationmaximization algorithm[END_REF] to the T1C image, which consistently has the highest resolution among the different MR channels.

The validation of the experiment consists in comparing real MRI with synthetic MRI generated solely based on a segmentation of the brain and the tumor. Quantitative assessment include Peak Signal to Noise Ratio (PNSR) to quantify the quality of the synthesis, Mean-Absolute-Error (MAE) to quantify bias, and finally mean Structural Similarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] (with default values) as an additional image reconstruction criterion. All these measures are classically used for assessment of image reconstruction or synthesis [START_REF] Cardoso | Template-Based Multimodal Joint Generative Model of Brain Data[END_REF], [START_REF] Manjón | Non-local MRI upsampling[END_REF]- [START_REF] Gómez | Joint Reconstruction of Multi-Contrast MRI for Multiple Sclerosis Lesion Segmentation[END_REF]. Whenever an average is computed, only the voxels which belong to a region-of-interest (ROI) are considered: the ROI is the whole brain mask as defined by the input label map (segmentation of healthy and pathological tissues).

Our image synthesis model is compared to the image synthesis procedure described in [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF]. Supplementary material, available in the supplementary files /multimedia tab, includes the user-specific coefficients chosen in our implementation of this procedure [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF].

For a fair comparison, synthetic MRI generated with [28] i) use the same input categorical label map as our method, and ii) are standardized as described in Section III-A, so that the average image intensity and the variance corresponds to the average and variance of the real image. Supplementary material, available in the supplementary files /multimedia tab, illustrates results obtained with [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF] based on tumor growth model probabilistic outputs instead. However, no quantitative evaluation can be performed on these results due to a lack of ground truth.

2) Results: A visual inspection of image synthesis results can be performed in axial views for a high-grade glioma in Figure 6. Axial, coronal and sagittal views are also presented in smaller displays in supplementary material, available in the supplementary files /multimedia tab. The model does not replicate artifacts or croppings present on some real cases. This is expected, since the image synthesis model is based on a label map, which was created from a high-resolution T1C image of a whole brain.

Quantitative results are reported for each MR channel in Figure 7 for PSNR, and in supplementary material, available in the supplementary files /multimedia tab, for other assessment measures. A summary of assessment measures is presented in Table II. PSNR, MAE, and mean SSIM are always higher for the proposed image synthesis method. Moreover, with the linearly increasing feedback weight, most results are improved with additional iterations. To confirm this statement, paired two-sample t-tests are performed independently for each MR channel. With 1% significance level, EMP with 2 iterations performs significantly better than the state-of-the-art [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF] for every criterion (PSNR/SSIM/MAE) and for every MR channel (T1/T1C/T2/FLAIR), which is not the case when using EMP with a single iteration. Moreover, with 1% significance level, from iteration 1 to iteration 3, EMP improves significantly the image between each iteration.

The influence of the feedback weight α t was studied. The Fig. 6: Qualitative evaluation of our synthesis method after 3 iterations on a high-grade case (HG-0009). Axial views. Top to bottom: synthetic MRI obtained with our method, ground truth (real MRI), synthetic MRI obtained with the method described in [START_REF]Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF] with the same categorical label map as input. Ground truth FLAIR image is partially dark, corresponding to the brain mask for which the FLAIR image is unknown, most likely for a faster MR acquisition. Interestingly, this does not happen on T1C images, and hence does not affect the label map for healthy brain tissues. Synthetic images are not affected since they rely only on the label map. Supplementary material, available in the supplementary files /multimedia tab, presents maps of absolute differences between synthesized images and ground truth. protocol is the following: different values are chosen for feedback weights, then Extended Modality Propagation is run for 10 iterations with α 1 = 0 and subsequent α t fixed at the chosen value. The highest PSNR (18.6) and highest mean-SSIM (0.63) are obtained for the smallest non-zero feedback weight (1/9). The variance of the PSNR is lower when a nonzero feedback weight is chosen, however the median PNSR is nearly the same for all feedback weight values. MAE increases when the weight increases. In the end, any small non-zero feedback weight is reasonable for image synthesis.

The influence of the number of iterations of EMP was studied. A high constant feedback weight, fixed at 2/3, is chosen to emphasize the effect of the number of iterations. The highest PSNR (18.7) and highest mean SSIM (0.64) are obtained at the end of the second iteration. MAE decreases notably with the number of iterations. In the end, a small number of iterations should be preferred, typically 2 to 3 iterations. A similar trend is observed with a linearly increasing feedback weight α t . [START_REF] Swanson | Quantifying the role of angiogenesis in malignant progression of gliomas: in silico modeling integrates imaging and histology[END_REF]. It was also used in 3D to simulate the impact of anti-angiogenic therapy in an atlas of a healthy brain [START_REF] Saut | A multilayer grow-or-go model for gbm: effects of invasive cells and antiangiogenesis on growth[END_REF]. Supplementary material, available in the supplementary files /multimedia tab, details the model and its implementation. Synthetic MRI based on the PIHNA model are shown in Figure 8.

IV. DISCUSSION

In the experiments presented in this paper, the iterative nature of the proposed image synthesis method allows to reach higher PSNR and synthesize qualitatively more realistic images. The optimal number of iterations seems to be about 2 or 3. A progressive increase of the feedback weight with the number of iterations, as suggested in [START_REF] Ye | Modality Propagation: coherent synthesis of subject-specific scans with data-driven regularization[END_REF] leads to better synthesis results.

Automatically setting the value for the feedback weight would require to consider the weight as an additional random variable together with an hyper prior which depends on the amount of spatial correlation of J.

Several MR channels are simultaneously synthesized, which leads to a more constrained patch match (for iterations t > 1), and ultimately to more consistent synthetic MRI. However, using a higher number of MR sequences implies that patches lie in a space of higher dimension, which tends to homogenize patch distances. Indeed, for high dimensional data, every patch is more or less a close neighbor of others, and the notion of nearest-neighborhood becomes less relevant. This could result in fuzzier synthetic images. Conversely, less consistent, yet potentially more textured, results could arise when trying to synthesize a single image modality.

EMP relies on a marginalization over several training cases, which has an averaging effect. Therefore, if labeling errors are assumed to be incoherent in the training cases, an incorrect labeling should have little impact, theoretically proportional to 1 over the number of atlases. However, if the target label map L is close to one of the incorrect segmentations, then the impact would be greater in practice.

One limitation of the proposed method is the running time: with 20 atlases, one iteration to simultaneously synthesize the 4 MR channels (for the whole brain) requires about 2 days of computation (16 cores, 100 GB of RAM) for 1 mm isotropic MRI, and about 3 hours of computation (3 cores, 20 GB of RAM) for 2 mm isotropic MRI. Typically, sulci are less visible on 2 mm synthetic MRI. Atlas selection could result in a decrease of the number of atlases, which could ultimately lead to finer details on 2 mm synthetic MRI. However, this would require a strategy to pre-select relevant atlases based on some similarity criterion between label maps.

Regarding running time in the perspective of generating a large database of synthetic cases, a region-of-interest enclosing the tumor could be defined such that the rest of the synthetic image is supplied by a template. Indeed, tumor growth models are based on a healthy atlas in which a tumor seed is placed, and a T1 template is available for the healthy atlas [START_REF] Mazziotta | A probabilistic atlas and reference system for the human brain: International Consortium for Brain Mapping (ICBM)[END_REF]. In this case, image intensity does not have to be optimized for voxels outside the ROI, since these voxels are assumed not to present intensity abnormalities.

Without resorting to tumor growth models, by only performing the patch-based synthesis on pathological regions, the extended modality propagation could be applied to any brain MR image, provided that a realistic label map of tumor compartments exists for that image. The generation of those tumor label maps could be either borrowed and transported from real pathological cases in the BRATS database or could be generated procedurally through an additional label map model. This would allow to create a high number of complex and detailed cases, with different locations of the same tumor label maps in the brain. The impact of this perspective would have to be carefully considered since some tumor patterns might be more or less likely in different brain regions.

All experiments were conducted with a label map as single input. Due to the iterative process described in this paper, it would be straightforward to provide as input an MR channel on top of the label map. Typically, if a contrast-enhanced T1 MR image is available on top of the label map, we could synthesize MR channels such as FLAIR or T2. Additional inputs are merely additional constraints for the patch match.

V. SUMMARY AND CONCLUSIONS

We presented a generative model for the synthesis of multimodal medical images of pathological cases by using an iterative multi-atlas patch-based algorithm. The method was illustrated on MRI synthesis of glioma-bearing brains. A single label maps allows to synthesize realistic images, with clear improvements with regard to the state-of-the-art. The approach is flexible and can be applied i) either to label maps of real cases to synthesize artifact-free MRI, ii) or to label maps simulated by any tumor growth model. The standardization of images is made transparent and reversible, so that synthetic data can be used for any application with any other dataset, including intensity inhomogeneity or different levels of Rician noise. Moreover, due to the strategy used to solve for the maximum a posteriori, estimation of image synthesis uncertainty is straightforward.

Such a generative model would allow the availability of large, public and annotated datasets of synthetic medical images obtained with virtual patients, which could prove useful for the validation or benchmarking of image processing algorithms. Promising results have been obtained on real data after training machine learning algorithms on synthetic data [START_REF] Geremia | Brain tumor cell density estimation from multi-modal MR images based on a synthetic tumor growth model[END_REF], which could be further improved by using a more realistic synthetic dataset, generated with the proposed image synthesis model. Finally, image synthesis of pathological cases bridges the gap between outputs of tumor growth model and medical images, which are the most clinically relevant observations. The development of image synthesis could benefit to interpretability and personalization of tumor growth models.

Fig. 1 :

 1 Fig. 1: The generative model aims at synthesizing subjectspecific MRI, provided a label map of a pathological brain.

Algorithm 1 :Fig. 4 :

 14 Fig. 4: Illustration of image synthesis uncertainty estimation.Along with the model prediction J(x) (synthesis mean), the knowledge of the model uncertainty is illustrated with the standard-deviation σ(x) (square-root of the synthesis variance). Please note that only the FLAIR MR sequence is shown here, but J(x) consists in fact in four synthesized MR channels, which are taken into account by σ(x) shown here. Uncertainty is more important at the boundary between cerebrospinal fluid and brain tissues (ventricles, sulci) and in the tumor core.

Fig. 7 :

 7 Fig. 7: Comparison to ground truth MRI in terms of PSNR for each MR sequence. From left to right: state-of-the-art [28], iterations 1, 2 and 3 of our method.

Fig. 8 :

 8 Fig. 8: Illustration of the modularity of our image synthesis model. Synthetic 1 mm isotropic MRI, obtained after 3 iterations, based on a categorical label map generated by the PIHNA tumor growth 3D model.

TABLE I :

 I Symmetry between likelihoods for segmentation & synthesis. Observed variables are displayed with grey background.

TABLE II :

 II Average assessment measures for image synthesis, based on the 2013 BRATS Evaluation dataset.

	T1	PSNR	mean SSIM	MAE
	State-of-the-art [28]	21.4	0.70	80.7
	Ours (iteration 1)	22.5	0.76	65.6
	Ours (iteration 2)	22.7	0.81	63.0
	Ours (iteration 3)	22.8	0.82	62.6
	T1C	PSNR	mean SSIM	MAE
	State-of-the-art [28]	21.9	0.77	75.2
	Ours (iteration 1)	22.0	0.76	70.1
	Ours (iteration 2)	22.1	0.81	68.6
	Ours (iteration 3)	22.2	0.82	68.0
	T2	PSNR	mean SSIM	MAE
	State-of-the-art [28]	16.5	0.67	97.2
	Ours (iteration 1)	18.5	0.69	76.0
	Ours (iteration 2)	18.7	0.74	73.3
	Ours (iteration 3)	18.8	0.75	73.0
	FLAIR	PSNR	mean SSIM	MAE
	State-of-the-art [28]	18.7	0.63	95.4
	Ours (iteration 1)	19.3	0.66	78.2
	Ours (iteration 2)	19.8	0.73	73.6
	Ours (iteration 3)	19.9	0.75	73.1
	C. Illustration using an advanced tumor growth model
	In this section, we use an advanced tumor growth model
	to generate probabilistic label maps, based on which we
	create a hard label map. The proposed image synthesis model
	is then used to generate synthetic MRI of virtual patients.
	The tumor growth model is a multi-population cell model
	driven by angiogenesis and vascularization. It is a variation
	of the Proliferation Invasion Hypoxia Necrosis Angiogenesis
	(PIHNA) model described in 1D in

Freely available online at http://www.nitrc.org/projects/tumorsim

ACKNOWLEDGEMENTS

Part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA 2011-291080 (on Biophysical Modeling and Analysis of Dynamic Medical Images).