Maïa Proisy

Bertrand Bruneau

Laurent Riffaud
email: laurent.riffaud@chu-rennes.fr

Pr L Riffaud

How ultrasonography can contribute diagnosis of craniosynostosis Comment l'échographie peut aider au diagnostic de craniosténose

Keywords: Craniosynostosis, suture, infant, ultrasound, skull

Introduction. The aim of this article was to provide an overview of ultrasound (US) techniques for the investigation of cranial sutures in infants.

Material and methods. We first describe a high-resolution sonography technique and its limitations. We then analyze the reliability, effectiveness and role of ultrasonography in routine practice using a PubMed literature review.

Results. Ten studies reported excellent correlations between ultrasonography and 3D-CT.

Cranial US for the diagnosis of a closed suture had 100% sensitivity in 8 studies and 86-100% specificity before the age of 12 months. Negative findings mean imaging investigation can be stopped. If ultrasonography confirms diagnosis, neurosurgical consultation is required. Thus, 3D-CT can be postponed until appropriate before surgery. Conclusion. Cranial suture ultrasound is an effective and reliable technique for the diagnosis of craniosynostosis. It has many advantages: it is fast and non-irradiating, and no sedation is required. It should be used as first-line imaging in infants below the age of 8-12 months when craniosynostosis is clinically suspected.

Introduction

Several studies have shown the effectiveness of ultrasound (US) in craniosynostosis, but it is not used for routine screening in many centers, as this depends on the preferences of referring physicians. The present study first discusses US techniques for the exploration of cranial sutures; we then review the literature, and finally propose a decision-tree for suspected craniosynostosis in infants.

Ultrasonography examination of cranial sutures

The sonographic aspect of normal major cranial sutures, in neonates and infants, was first described by Soboleski et al. [START_REF] Soboleski | Sonography of normal cranial sutures[END_REF]. High-resolution sonography of abnormal cranial sutures was described 1 year later by the same author [START_REF] Soboleski | Highresolution sonography of the abnormal cranial suture[END_REF]. A high frequency linear transducer with a minimum 12.5 MHz (and ideally >17.5 MHz) should be used, with parameters adjusted depending on the device. The sonographic examination may begin at the anterior fontanelle, with a coronal orientation. A transverse orientation, perpendicular to the plane of the long axis of the suture, has to be maintained along the length of the suture. Coronal, sagittal, lambdoid and metopic sutures must be followed along their whole length.

A normal suture is defined as a hypoechoic gap between two hyperechoic bony plates. The sagittal suture usually has an end-to-end appearance; the lambdoid and coronal sutures usually have a beveled or overlapped appearance [START_REF] Soboleski | Sonography of normal cranial sutures[END_REF] (Fig. 1). There was no correlation between suture width or thickness and patient age between 0 and 5 months in Soboleski et al.'s study [START_REF] Soboleski | Sonography of normal cranial sutures[END_REF]. A synostosed suture shows loss of the echo-poor fibrous gap between bony plates [START_REF] Soboleski | Highresolution sonography of the abnormal cranial suture[END_REF]. An osseous ridge may be observed, in addition to the loss of the normal hypoechoic gap (Fig. 2). Closure may be partial or complete (Fig. 3).

J o u r n a l P r e -p r o o f

Ultrasound imaging is often described as an "operator-dependent" technique. Although education of radiologists and a short training module is recommended, cranial suture US is relatively easy to perform. Furthermore, it has been shown that there is no interobserver variability with this method [START_REF] Regelsberger | Highfrequency ultrasound confirmation of positional plagiocephaly[END_REF]. Analysis of coronal sutures may be a little more difficult because of the overlap and the difficulty of maintaining a transverse orientation to the suture while scanning.

Validity of the technique

A comprehensive literature review of the National Library of Medicine (PubMed) database was performed to identify relevant studies. It was conducted in February 2019 and analyzed articles published up to that date, with no starting date limit. We used the following key-words: "craniosynostosis diagnosis ultrasound sensitivity" and "cranial suture ultrasound diagnosis craniosynostosis". Studies were selected if their Abstracts were written in English and reported an original article comparing US to a gold standard examination (CT and/or X-ray and/or clinical follow-up) for exploring cranial suture patency. Ultrasound sensitivity and specificity were required to be reported in the study.

The effectiveness and reliability of US examination in children with suspected craniosynostosis, in comparison with radiography, CT or clinical follow up, has been demonstrated in several studies [START_REF] Sze | Ultrasound screening of the lambdoid suture in the child with posterior plagiocephaly[END_REF][START_REF] Regelsberger | Ultrasound in the diagnosis of craniosynostosis[END_REF][START_REF] Simanovsky | Effectiveness of ultrasonographic evaluation of the cranial sutures in children with suspected craniosynostosis[END_REF][START_REF] Krimmel | Value of highresolution ultrasound in the differential diagnosis of scaphocephaly and occipital plagiocephaly[END_REF][START_REF] Alizadeh | Diagnostic accuracy of ultrasonic examination in suspected craniosynostosis among infants[END_REF][START_REF] Linz | Occipital plagiocephaly: unilateral lambdoid synostosis versus positional plagiocephaly[END_REF][START_REF] Rozovsky | Cranial Ultrasound as a First-Line Imaging Examination for Craniosynostosis[END_REF][START_REF] Pogliani | Cranial ultrasound is a reliable first step imaging in children with suspected craniosynostosis[END_REF][START_REF] Hall | Effectiveness of screening for craniosynostosis with ultrasound: a retrospective review[END_REF][START_REF] Proisy | Ultrasonography for the diagnosis of craniosynostosis[END_REF], showing high sensitivity and specificity (Table 1). In addition, the capacity to determine complete or incomplete closure of the sutures also reaches 100% [START_REF] Proisy | Ultrasonography for the diagnosis of craniosynostosis[END_REF]. Inconclusive US examination is very infrequent: it may be due to poor cooperation in some studies [START_REF] Rozovsky | Cranial Ultrasound as a First-Line Imaging Examination for Craniosynostosis[END_REF] or to relatively advanced age (12-18 months) [START_REF] Pogliani | Cranial ultrasound is a reliable first step imaging in children with suspected craniosynostosis[END_REF][START_REF] Hall | Effectiveness of screening for craniosynostosis with ultrasound: a retrospective review[END_REF]. The less conclusive results in Krimmel's study [START_REF] Krimmel | Value of highresolution ultrasound in the differential diagnosis of scaphocephaly and occipital plagiocephaly[END_REF] can be explained by the inclusion of inconclusive US findings as false-positive and false-negative results in their statistics. Alizadeh et al. [START_REF] Alizadeh | Diagnostic accuracy of ultrasonic examination in suspected craniosynostosis among infants[END_REF] had one false J o u r n a l P r e -p r o o f negative case in a 7 month-old boy, but did not provide an explanation for this. Pogliani et al. [START_REF] Pogliani | Cranial ultrasound is a reliable first step imaging in children with suspected craniosynostosis[END_REF] reported 2 false positives during study start-up, and 1 false negative case of bi-temporal suture closure missed by US: they did not evaluate temporal sutures either for technical difficulties or for the lack of definitive literature data about their significance in true craniosynostosis.

Role of US in diagnosis strategy (Fig.4)

Diagnosis of positional plagiocephaly can be assessed on physical examination, based on typical features [START_REF] Linz | Occipital plagiocephaly: unilateral lambdoid synostosis versus positional plagiocephaly[END_REF] such as trapezoidal head shape, mastoid bulge, skull base tilt and pushedforward ear position. However, clinical differentiation between isolated unilateral lambdoid and coronal synostosis can be difficult, especially for general practitioners or pediatricians. Rozovsky et al. [START_REF] Rozovsky | Cranial Ultrasound as a First-Line Imaging Examination for Craniosynostosis[END_REF] recommended an algorithm for diagnosis or exclusion of non-syndromic craniosynostosis by pediatricians or general practitioners: when craniosynostosis is clinically suspected, a cranial US examination should be performed first. If negative, there is no need for further imaging studies. If cranial US confirms diagnosis or is indecisive, then neurosurgical consultation is required.

The main advantages of US examination are that it is a fast non-irradiating technique which does not require sedation. It is easy to carry out, facilitating early diagnosis. With US as first-line imaging procedure, 3D-CT may be performed subsequently in case of uncertain US clinical situations. In case of positive US findings in preoperative planning, 3D-CT may also be postponed before surgery at the surgeon's discretion. Considering the potential side-effects of ionizing radiation, sedation risks and costs, some authors even suggest that operative correction of isolated sagittal synostosis can proceed without any radiological investigations, unless clinical examination exhibits atypical features [START_REF] Fearon | The Diagnosis and Treatment of Single-Sutural Synostoses: Are Computed Tomographic Scans Necessary?[END_REF][START_REF] Agrawal | Diagnosis of isolated sagittal synostosis: are radiographic studies necessary? Child's[END_REF]. In our opinion, diagnosis still needs confirming, J o u r n a l P r e -p r o o f and cranial US can be a good alternative to 3D-CT. Clinical examination may occasionally be insufficient even for neurosurgeons. For example, some patients have a non-synostotic elongated cranium that suggests scaphocephaly [START_REF] Proisy | Ultrasonography for the diagnosis of craniosynostosis[END_REF][START_REF] Baumgartner | Nonsynostotic scaphocephaly: the so-called sticky sagittal suture[END_REF]; these children should not be operated on and, if they are treated early by molding helmet therapy (at < 12 months of age), head shape may become normalized. Also, breech babies are known to have a mild skull deformation (dolichocephaly) known as "breech head", which may mimic scaphocephaly [START_REF] Kasby | The breech head and its ultrasound significance[END_REF]. Finally, and regarding these last two points, when craniosynostosis is clinically suspected we consider it necessary to confirm diagnosis systematically on imaging.

According to our experience and the literature, US is limited by age, with the upper limit around 8-12 months [START_REF] Regelsberger | Highfrequency ultrasound confirmation of positional plagiocephaly[END_REF]. Reliability decreases with age because of hair growth, skull growth and thinner sutures. As most children with cranial deformity are usually seen in the first month of life, we recommend carrying out US examination as soon as possible when diagnosis is suspected, and ideally before the age of 8 months. Finally, and given that physiological closure occurs much earlier in the metopic suture than in the others (normally between the age of 3 and 9 months [START_REF] Vu | The timing of physiologic closure of the metopic suture: a review of 159 patients using reconstructed 3D CT scans of the craniofacial region[END_REF]), complete fusion of the metopic suture does not necessary mean trigonocephaly. Typical triangular frontal deformation and hypotelorism should also be associated with complete fusion of the metopic suture for trigonocephaly to be diagnosed with certainty.

Conclusion

US examination of cranial sutures has been shown in the literature to be an effective and reliable technique. Cranial US has many advantages: it is fast and easy to implement and nonirradiating, and infants do not need to undergo sedation. It may serve as first-line imaging when a clinical diagnosis of craniosynostosis is suspected. It allows the 3D-CT examination, if required, to be postponed to an appropriate time before surgery. Given the virtual impossibility of

Figure 1 :Figure 2 :

 12 Figure 1: Ultrasonographic appearance of normal sutures. A = appearance of normal sagittal

Figure 3 :

 3 Figure 3: Ultrasonographic appearance of sutures (A) and lateral 3D-CT view (B) in a 3 month-

Figure 4 :

 4 Figure 4: Decision-tree for diagnosis of craniosynostosis.

This work did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

J o u r n a l P r e -p r o o f assessing suture patency after a certain age (around 8-12 months), the earlier the US scan is performed the more reliable the examination.

Disclosure of interest

The authors declare that they have no competing interests.


Jaakko Nevalainen

Laurent Guyon

Firstly I would like to thank my supervisors

Chapter 1 presents an overview of undirected graphical models in general and Markov models in particular. Chapter 2 presents inference for variables in hidden Markov random fields (MRFs) while Chapter 3 presents inference for the parameters of Gaussian MRFs. Chapter 4 outlines the four applications and how the Markov model framework is utilised in each case. For each application, the associated publication is also provided.

In Publication 1, hidden Markov models (HMMs) are used to achieve an alignment and classification of time series data. Publications 2 and 3 concern inference with hidden MRFs to obtain a segmentation of both digital image data and network data respectively. Spatial analysis with Gaussian MRFs is presented in Publication 4. We show that our particular use of Markov models in each of our applications enables us to achieve our aims.

Tiivistelmä

Tässä väitöskirjassa esitellään neljä Markovin mallien bioinformatiikan sovellusta. Väitöskirjassa laajennetaan Markovin mallien käyttöä monelta eri sovellusalalta kerätyn biologisen datan matemaattiseen ja tilastolliseen analyysiin. Malleja hyödynnetään genomitasolla aikasarjoissa ja verkkodatassa sekä solutasolla soluviljelyn ja in vivo -kudoksen mikroskooppikuva-aineistossa.

Kiinnostuksen kohteina tässä tutkimuksessa ovat geeneistä, soluista ja pikseleistä koostuvat objektien joukot. Työssä käytetään hyväksi näiden joukkojen sisällä vallitsevia spatiaalisia, temporaalisia tai funktionaalisia assosiaatiota olettaen, että toisiaan lähempänä olevat objektit ovat keskenään vahvemmin yhteyksissä kuin toisistaan kauempana olevat. Tämä oletus mahdollistaa tehokkaan päättelyn Markov-malliperheessä sekä ehdollisten riippumattomuuksien esittämisen satunnaismuuttujien välillä, joita kuvataan suuntaamattomissa graafeissa pisteillä ja viivoilla.

Luvussa 1 tarkastellaan yleisellä tasolla suuntaamattomia malleja, joista tarkemmin esitellään Markovin mallit. Luvussa 2 käsitellään muuttujia koskevaa päättelyä piilo-Markovin satunnaiskentissä (MRF) ja luvussa 3 keskitytään Gaussisten MRF -mallien parametreja koskevaan päättelyyn. Luvussa 4 esitellään neljä tutkimuksen sovellusta ja miten Markovin malleja on niissä laajennettu ja hyödynnetty. Sovelluksiin liittyvät osatyöt muodostavat väitöskirjan viimeisen luvun.

Osatyössä 1 Markovin piilomalleja (HMM) käytetään aikasarjadatan kohdennukseen ja luokitteluun. Osatöissä 2 ja 3 käsitellään päättelyä piilo-MRF-malleissa niiden hyödyntämiseksi digitaalisen kuva-ja verkkodatan segmentoimisessa. Osatyössä 4 käsitellään Gaussisten MRF -mallien käyttöä spatiaalisessa analyysissa. Huolellisesti rakennettujen Markovin mallien käyttö osoittautuu kaikissa sovelluksissa keskeiseksi osatutkimusten tavoitteiden saavuttamiseksi.

Résumé

Dans cette thèse nous présentons quatre applications en bioinformatique avec des modèles de Markov. Plus précisément, nous étendons ces modèles à l'analyse statistique et mathématique de données biologiques. Les données que nous étudions viennent de différentes sources. Nous considérons les applications au niveau génomique avec des séries temporelles et des données de réseau ainsi que les applications au niveau cellulaire avec des données d'images de microscopie de cultures cellulaires et de tissus in vivo.

Les ensembles d'objets tels que gènes, cellules ou pixels, présentent un intérêt dans leur intégralité. Nous utilisons des associations spatiales, temporelles ou fonctionnelles au sein de ces ensembles et nous supposons que les objets qui sont plus proches les uns des autres sont plus fortement liés que ceux qui sont plus éloignés. Cela permet une inférence efficace dans le cadre des modèles de Markov et est codé par des indépendances conditionnelles entre variables qui sont représentées par des sommets et arêtes dans un graphe non orienté.

Le chapitre 1 présente une vue d'ensemble des modèles graphiques non orientés en général et des modèles de Markov en particulier. Le chapitre 2 présente l'inférence des variables des champs aléatoires de Markov (MRFs) cachés tandis que le chapitre 3 présente l'inférence des paramètres des MRFs de Gauss. Le chapitre 4 expose les quatre applications traitées et comment le cadre des modèles de Markov est utilisé dans chaque cas. Pour chacune des applications, la publication associée est fournie.

Introduction

Markov models are a type of undirected graphical model, a collection of random variables whose conditional independence properties can be represented by an undirected graph. They are used for many different applications in a broad range of areas including the analysis of time series data, image data, network data and spatial data among many others [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF][START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF], Blake et al. 2011[START_REF] Banerjee | Hierarchical Modeling and Analysis for Spatial Data[END_REF]. Markov (1856Markov ( -1922)), such models are especially popular since the Markov structure allows for complex conditional independences to be modelled while still allowing for efficient inference. Moreover, the graphical representation of the conditional independence properties provides a valuable and intuitive visualisation of the model.

Named after Russian mathematician Andrei

Below we define undirected graphical models in general before defining Markov models in particular. We then consider the Hammersley-Clifford Theorem, a fundamental result allowing for the factorisation of the model density. The content of this chapter is primarily based on the presentations of [START_REF] Whittaker | Graphical Models in Applied Multivariate Statistics[END_REF] and [START_REF] Lauritzen | Graphical Models[END_REF].

Graphical Models

We use p to denote a probability density or probability mass function and use the argument to identify the random variable concerned. Hence for the random variables X and Z we write p(x) to represent the marginal density function, p(x, z) to represent the joint density function and p(x|z) to represent the conditional density function. A random variable X takes values in its support Ω X .

For collections of random variables indexed by scalar subscripts, {X 1 , X 2 , . . .} we use set notation in the subscript, X A = {X i | i ∈ A} along with standard set notation operations.

for all x ∈ Ω X , y ∈ Ω Y and z ∈ Ω Z such that p(z) > 0. Definition 1.2 Undirected Graph. † An undirected graph G is a set of vertices V and a set of edges E, written G = (V, E) where E is a subset of the set V × V of pairs of distinct vertices such that (i, j) ∈ E ⇔ (j, i) ∈ E for all i, j ∈ V. Definition 1.3 Conditional Independence Graph. † The conditional independence graph of the collection of random variables X V is the undirected graph G = (V, E) where (i, j) / ∈ E ⇔ X i ⊥ ⊥ X j | X V\{i,j} (1.1)
for all i = j ∈ V.

We have defined a conditional independence graph in terms of the pairwise Markov property (1.1). Given a collection of random variables with a set of conditional independence properties, we are able to construct the corresponding conditional independence graph. Under certain conditions that are given below, we are able to deduce further conditional independence properties of the model directly from the corresponding graph.

Definition 1.4 Paths and Separation. † In an undirected graph G, a path between two vertices i, j ∈ V is a sequence a 1 , a 2 , . . . , a K such that a 1 = i, a K = j and (a k-1 , a k ) ∈ E for k = 2, 3, . . . , K. Two vertices i, j ∈ V are separated by the subset C ⊆ V\{i, j} if and only if all paths between i and j contain at least one member of C.

Theorem 1.1 Separation Theorem. Let the collection of random variables X V have density function p(x V ) > 0 for all x V ∈ Ω X V and conditional independence graph G. For disjoint subsets of vertices A, B, C ⊂ V, if in the conditional independence graph G each vertex in A is separated from each vertex in B by the subset C, then

X A ⊥ ⊥ X B | X C . (1.2)
Proof of Theorem 1.1. A proof is given by [START_REF] Lauritzen | Graphical Models[END_REF]. † These definitions are taken directly from my Master of Philosophy thesis [START_REF] Robinson | Alignment of Time Course Microarray Data with Hidden Markov Models[END_REF]).

(i) = {j ∈ V | (i, j) ∈ E}. A clique is a subset of vertices, C ⊂ V such
that for all i ∈ C, if j ∈ C\i then j ∈ ν(i). A clique C is maximal if there exists no other clique

C * ⊆ V such that C ⊂ C * .
Each vertex is a clique and each pair of neighbouring vertices is also a clique. If a pair of neighbouring vertices have a mutual neighbour, then these three vertices are another clique. If all members of a clique don't have a mutual neighbour not already in the clique, then the clique is maximal.

Definition 1.6 Markov Model. The collection of random variables X V with conditional independence graph G is a Markov model if

X i ⊥ ⊥ X V\{i,ν(i)} | X ν(i) (1.3)
for all i ∈ V. Trees. An undirected graph G is a tree if there exists a unique path between each pair of vertices i, j ∈ V such that no vertex is contained in the path more than once. Markov models. Note that the vertices are labelled by the corresponding random variables. The major differences between the models in Figure 1.1 are whether they are trees or not (above and below) and whether the Markov components are observed or 'hidden' (left and right). In the 'hidden' models it is taken that the random variables Z V are observed while the random variables X V , constituting the Markov component of the model, are unobserved or missing. The defining difference between a Markov chain or hidden Markov model (HMM) as against a Markov random field (MRF) is that the latter is not a tree.

These differences determine the possible inference that can be achieved in the different models.

It is the tree structure of both Markov chains and HMMs that allows for relatively straightforward inference of both the hidden variables and model parameters [START_REF] Robinson | Alignment of Time Course Microarray Data with Hidden Markov Models[END_REF]). In the following we consider inference for the variables in hidden MRFs (Chapter 2) and the model parameters of Gaussian MRFs (Chapter 3). Note that in all of these cases the conditional independence graph itself is assumed to be known. Furthermore the graphs are also assumed to be sparse, that is, there are few edges between the vertices compared to all those possible. This sparsity is implicit in the notion of a local neighbourhood and the local Markov property.

Markov chain

X a X b X c Markov random field (MRF) X g X h X i X d X e X f X a X b X c Hidden Markov model (HMM) X a X b X c Z a Z b Z c Hidden MRF X g X h X i X d X e X f X a X b X c Z g Z h Z i Z d Z e Z f Z a Z b Z c
(x V ) > 0 for all x V ∈ Ω X V then p(x V ) ∝ C∈M F C (x C )
where M is the set of all maximal cliques in the graph and the factors F C (x C ) are strictly positive.

Proof of Theorem 1.4. Details of a proof are given by [START_REF] Lauritzen | Graphical Models[END_REF].

The factors F C (x C ) must be strictly positive but are otherwise arbitrary functions of their arguments. The Hammersley-Clifford Theorem is named after the unpublished paper by [START_REF] Hammersley | Markov Fields on Finite Graphs and Lattices[END_REF].

Theorem 1.5 The Hammersley-Clifford Theorem and the Markov Properties. If the collection of random variables X V has density function p(x V ) > 0 for all x V ∈ Ω X V then all three Markov properties and the Hammersley-Clifford Theorem are equivalent.

Proof of Theorem 1.5. A proof is given by [START_REF] Lauritzen | Graphical Models[END_REF].

This equivalency and the fact that the Hammersley-Clifford Theorem was not originally published means that it appears in many different forms [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems (with Discussion)[END_REF][START_REF] Lauritzen | Graphical Models[END_REF][START_REF] Cressie | Statistics for Spatio-Temporal Data[END_REF]. We follow [START_REF] Lauritzen | Graphical Models[END_REF] as this is the most convenient for us. It is the three Markov properties and the Hammersley-Clifford Theorem that allow for efficient inference so that Markov models can be gainfully used in practice.

Chapter 2

Inference for Hidden MRFs

Markov random fields (MRFs) have been extensively utilised in computer vision and in particular for digital image segmentation (Blake et al. 2011). Segmentation is the task of assigning a label to each pixel so that the image is partitioned into regions of pixels corresponding to mutually relevant features. Many different labelling problems can be approached as an inference problem in a suitably defined hidden MRF.

We consider inference for the random variables X V in a hidden MRF. In general such inference in graphical models is NP-hard [START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF]. However under certain conditions exact inference is computationally feasible for hidden MRFs while in other cases there exist efficient algorithms for approximate inference.

We first review the binary labelling problem as an inference problem for hidden MRFs. We present the energy minimisation framework and the use of graph cuts to find a computationally efficient and exact solution to the corresponding inference problem. Then we consider the multilabel problem and the established α-expansion algorithm for approximate energy minimisation.

Energy Function

For each i ∈ V, let X i be the random variable for the unobserved label we aim to infer and let Z i be the random variable for the observed data. Let the collection of these random variables {X V , Z V } be an hidden MRF with a conditional independence graph G, an example of which is given in the bottom right of Figure 1.1. In this case, the model has a regular grid structure but this is not necessary in general, just that each Z i has only X i as a unique neighbour. We consider that the random variables are simply indexed by scalar values for generality and convenience.

CHAPTER 2. INFERENCE FOR HIDDEN MRFS

The labelling task set as an inference problem is, given the hidden MRF {X V , Z V } and having observed the data Z V = z V , find the maximum a posteriori labels

xV = argmax x V p(x V |z V ) = argmax x V p(x V , z V ) p(z V ) = argmax x V p(x V , z V ).
For notational convenience let G be the conditional independence graph of only the random variables X V . That is, following our example graph in the bottom right of Figure 1.

1, G = (V, E)
is the graph given in the bottom left of Figure 1.1. Recall we assume that all Markov models we consider have a strictly positive density function and hence we use the Hammersley-Clifford Theorem (Theorem 1.4) to write

p(x V , z V ) ∝ i∈V F i (x i , z i ) (i,j)∈E F (i,j) (x i , x j )
where F i (x i , z i ) and F (i,j) (x i , x j ) are arbitrary, strictly positive functions of their arguments.

Definition 2.1 Energy Function. Given the hidden MRF {X V , Z V } with density function

p(x V , z V ) > 0 for all x V ∈ Ω X V and z V ∈ Ω Z V , and conditional independence graph G = (V, E)
only corresponding to the random variables X V , the associated energy function is

E(x V ) = i∈V φ i (x i ) + (i,j)∈E φ (i,j) (x i , x j ) (2.1)
with energy potentials

φ i (x i ) = -log(F i (x i , z i )) and φ (i,j) (x i , x j ) = -log(F (i,j) (x i , x j )).
This is just rewriting the factorisation given by Theorem 1.4 for a hidden MRF although note that we have not explicitly included z V in the energy function (2.1) as these terms are now accounted for in the subscripts of the energy potentials. Now we have that p(x V , z V ) ∝ exp{-E(x V )} and so the maximum a posteriori /minimum energy labels are therefore

xV = argmax x V p(x V |z V ) = argmin x V E(x V ).
The use of an energy function in computer vision goes back to [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], who established the connection with Gibbs distributions (Josiah Willard Gibbs (1839-1903) was an American mathematician and physicist). [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF] also coined the term 'Gibbs sampling' in their approach to energy minimisation, which was superseded by the iterated conditional modes (ICM) algorithm in computer vision [START_REF] Besag | On the Statistical Analysis of Dirty Pictures (with Discussion)[END_REF]).

The current use of energy minimisation in computer vision was brought about by the development of techniques based on graph cuts in the early 2000s, which allowed for much greater computational efficiency and improved performance as against ICM [START_REF] Szeliski | A Comparative Study of Energy Minimization Methods for Markov Random Fields with Smoothness-Based Priors[END_REF][START_REF] Kappes | A Comparative Study of Modern Inference Techniques for Structured Discrete Energy Minimization Problems[END_REF]. Note that in the hidden MRF given in Figure 1.1 the maximal cliques are pairs of vertices and for arbitrary graphs the maximal cliques may be larger than just pairs. However when the graph is still sparse a general approach is to just consider pairwise cliques in order to have an approximate energy function that can be optimised in practice [START_REF] Boykov | Fast Approximate Energy Minimization via Graph Cuts[END_REF]).

Graph Cuts

We first consider the binary label minimisation problem of (2.1), that is, where x i ∈ {0, 1} for all i ∈ V. Rewrite the energy potentials

φ i (x i ) =      θ i;0 if x i = 0 θ i;1 if x i = 1 and φ (i,j) (x i , x j ) =                  θ ij;00 if x i = x j = 0 θ ij;01 if x i = 0, x j = 1 θ ij;10 if x i = 1, x j = 0 θ ij;11 if x i = x j = 1
for all i ∈ V and (i, j) ∈ E. Hence we can rewrite the energy function (2.1) as

E(x V ) = i∈V θ i;1 x i + θ i;0 (1 -x i ) + (i,j)∈E θ ij;11 x i x j + θ ij;10 x i (1 -x j ) + θ ij;01 (1 -x i )x j + θ ij;00 (1 -x i )(1 -x j ) . (2.
2)

The approach used to minimise the energy (2.2) is to recast the problem as one of finding the maximum flow in a suitably defined network. Then, provided the energy potentials satisfy certain conditions, the minimum energy solution exists and can be found using graph cuts. The following material on network flow is based on the presentation of [START_REF] Bondy | Graph Theory[END_REF].

Definition 2.2 Directed Graph. A directed graph - → G is a set of vertices V and a set of edges - → E , written - → G = (V, - → E ), where - → E is a subset of the set V × V of pairs of distinct vertices.
A directed graph is an undirected graph (Definition 1.2) without the condition that if

(i, j) ∈ - → E then (j, i) ∈ - → E for all i, j ∈ V.
That is, it may be the case that only one of the two possible edges between vertices i and j are in the edge set -→ E . Hence there is a direction to each edge that allows for a different definition of neighbouring vertices.

→ ν (i) = {j ∈ V | (i, j) ∈ - → E } and ← - ν (i) = {j ∈ V | (j, i) ∈ - → E }. Definition 2.4 Network. A network N is a directed graph - → G with specified vertices s, t ∈ V,
known as the source and the sink respectively, and associated capacity function c :

V × V → [0, ∞) where c(i, j) = 0 if (i, j) / ∈ - → E , written N = (V, - → E , c).
A network is a particular type of directed graph where each edge has an associated capacity allowing for a flow between the specially designated source s and sink t vertices to be defined.

Definition 2.5 Flow in a Network. A flow in a network N is a function f : V × V → [0, ∞) such that 0 ≤ f (i, j) ≤ c(i, j) for all (i, j) ∈ - → E and j∈ - → ν (i) f (i, j) = j∈ ← -ν (i) f (j, i)
for all i ∈ V\{s, t}. The value of the flow is

val(f ) = j∈ - → ν (s) f (s, j) - j∈ ← -ν (s)
f (j, s).

(2.

3)

The flow on an edge is less than or equal to the capacity of that edge and for any network there exists a zero flow where f (i, j) = 0 for all (i, j) ∈ -→ E . Flow is conserved for all vertices except the source s and sink t. That is, the flow into any vertex i ∈ V\{s, t} (the sum of the flows on the edges (j, i) ∈ -→ E ) is equal to the flow out of i (the sum of the flows on the edges (i, j) ∈ -→ E ).

Lemma 2.1 Value of the Flow. For a flow f in a network N the value of the flow is equivalently

val(f ) = i∈ ← -ν (t) f (i, t) - i∈ - → ν (t)
f (t, i).

(2.4)

Proof of Lemma 2.1. A proof is given by [START_REF] Bondy | Graph Theory[END_REF].

That is, the value of the flow is defined to be the flow out of the source minus the flow into the source (2.3), which is equal to the flow into the sink minus the flow out of the sink (2.4).

Definition 2.6 Graph Cut. For a subset of vertices S ⊂ V in a network N with s ∈ S and t / ∈ S, a cut is the set of edges

[S] = {(i, j) ∈ - → E | i ∈ S, j / ∈ S}
and the value of the cut is

val[S] = (i,j)∈[S]
c(i, j).

A cut in a network N is the set of edges from the subset S to its complement that effectively 'cut' the network in two and hence 'graph cuts'. Proof of Theorem 2.2. A proof is given by [START_REF] Bondy | Graph Theory[END_REF].

Given a network N , the problem is to find a flow f and a cut [S] such that val(f

) = val[S].
That is, the value of the flow is maximised and the value of the cut is minimised. There exist classic algorithms for this task [START_REF] Ford | Maximal Flow Through a Network[END_REF]).

We will now consider how the maximum flow/minimum cut framework can be applied to the energy minimisation problem with a binary label energy function (2.2). First, we create a suitable network using the conditional independence graph of X V .

Algorithm 2.1 Network Creation Algorithm. Given the collection of random variables

X V with conditional independence graph G = (V, E), create a corresponding network N * = (V * , - → E * , c * )
by defining:

1. The set of vertices V * = V ∪ {s, t} where s is the source and t is the sink.

The edge set

- → E * = (s, i), (i, s), (t, i), (i, t) | i ∈ V ∪ (i, j), (j, i) | (i, j) ∈ E .
3. The capacity function

c * (i, j) = ω ij ≥ 0 for all (i, j) ∈ - → E * with c * (i, s) = 0 for all i ∈ V and c * (t, j) = 0 for all j ∈ V.
Consider that any cut [S] in the network N * can be viewed as a binary labelling

x i =      0 if i ∈ S 1 if i / ∈ S (2.5)
for all i ∈ V and consider the energy function

E * (x V ) = i∈V ω si x i + ω it (1 -x i ) + (i,j)∈E ω ij (1 -x i )x j + ω ji x i (1 -x j ) .
(2.6) Theorem 2.3 Graph Cut Energy. Given the collection of random variables X V with conditional independence graph G = (V, E), the value of a cut [S] in the corresponding network N * is given by (2.6) when utilising the labelling (2.5). That is, val[S] = E * (x V ).

Proof of Theorem 2.3. for all (i, j) ∈ E.

val[S] = (i,j)∈[S] ω ij = (i,j)∈ - → E * ω ij I{i ∈ S, j / ∈ S} = (s,i)∈ - → E * ω si I{i / ∈ S} + (i,t)∈ - → E * ω it I{i ∈ S} + (i,j)∈ - → E * ω ij I{i ∈ S\s, j / ∈ S ∪ t} = i∈V ω si I{i / ∈ S} + ω it I{i ∈ S} + (i,j)∈E ω ij I{i ∈ S, j / ∈ S} + ω ji I{i / ∈ S,
For the binary label energy function (2.2), constants can be added and subtracted to the potentials so that the energy function remains unchanged, that is the function still gives the same output for the same input. This is known as a reparameterisation of the energy function and if the potentials are submodular, then the energy function (2.2) can be reparameterised so that the potentials are in the form (2.7) in a finite series of basic prescribed operations [START_REF] Kolmogorov | What Energy Functions Can Be Minimized via Graph Cuts?[END_REF]. There are standard energy minimisation algorithms based on graph cuts for the binary label, submodular potential problem [START_REF] Boykov | An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision[END_REF]). Below, we consider the case where there are more than 2 labels.

Alpha-Expansion Algorithm

Now consider the case where we have more than 2 labels, that is x i ∈ L for all i ∈ V where L is a set of more than 2 elements. Consider the original energy function (2.1) and rewrite

φ i (x i = l) = θ i;l and φ (i,j) (x i = l, x j = k) = θ ij;lk
for all i ∈ V, (i, j) ∈ E and l, k ∈ L so that [START_REF] Boykov | Fast Approximate Energy Minimization via Graph Cuts[END_REF]).

E(x V ) = i∈V l∈L θ i;l I{x i = l} + (i,j)∈E l∈L k∈L θ ij;lk I{x i = l, x j = k}. ( 2
However, we can approximate the minimum energy solution using an algorithm that is based on graph cuts [START_REF] Boykov | Fast Approximate Energy Minimization via Graph Cuts[END_REF][START_REF] Boykov | An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision[END_REF][START_REF] Kolmogorov | What Energy Functions Can Be Minimized via Graph Cuts?[END_REF]. The main idea is to recast the multi-label problem in terms of a binary label problem by considering an initial labelling and then a binary choice for each variable: either change label or stay the same. In order to properly determine the problem, 'changing label' needs to be defined and one option is to only allow a change to a specified label, α ∈ L. The following material is based on the presentation of [START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF].

Let the initial labelling be x * V . We construct a binary label problem with dummy variables Y V such that y i ∈ {0, 1} for all i ∈ V. Consider that from the initial labelling x * V , the dummy variables Y V and a specified label α ∈ L, we obtain a new labelling xV given by

xi =      x * i if y i = 0 α if y i = 1
for all i ∈ V. Hence each variable either keeps its initial label x * i or changes to α depending on the value of y i . By substituting xV into (2.8) we obtain the equivalent energy function

Ẽ(y V ) = i∈V φi (y i ) + (i,j)∈E φ(i,j) (y i , y j ) (2.9) where φi (y i ) =      θ i;x * i if y i = 0 θ i;α if y i = 1 and φ(i,j) (y i , y j ) =                  θ ij;x * i x * j if y i = y j = 0 θ ij;x * i α if y i = 0, y j = 1 θ ij;αx * j if y i = 1, y j = 0 θ ij;αα if y i = y j = 1.
That is, we have a binary label energy function (2.9) for the binary label problem for Y V . By construction, Ẽ(y V ) = E(x V ) and hence finding the optimal labels ŷV is equivalent to finding the optimal labels xV in the restricted space of only allowing each variable to keep its initial label or change to α. This is a constrained energy minimisation problem where we are only (potentially)

increasing the number of vertices with the label α and hence known as 'alpha-expansion'. In practice, all of the possible labels are considered iteratively and multiple times. 2. Solve the binary label problem for Y V using graph cuts as outlined in Section 2.2.

3. Update the labelling and repeat for another label until some convergence criteria is satisfied.

Recall that in order to use graph cuts for the binary label problem, the energy potentials must be submodular. The potentials for the above binary label energy (2.9) are submodular when

θ ij;αα + θ ij;x * i x * j ≤ θ ij;x * i α + θ ij;αx * j
for all (i, j) ∈ E. This condition is satisfied when the energy potentials are a metric.

Definition 2.9 Semi-metric and Metric. A function

D : V × V → [0, ∞) is a semi-metric if D(i, j) = 0 ⇔ i = j and D(i, j) = D(j, i)
for all i, j ∈ V, and is a metric if in addition

D(i, k) ≤ D(i, j) + D(j, k)
for all i, j, k ∈ V.

However, even when the energy potentials are only a semi-metric, the α-expansion algorithm may still be gainfully used in practice [START_REF] Boykov | Fast Approximate Energy Minimization via Graph Cuts[END_REF].

Chapter 3

Inference for Gaussian MRFs

Gaussian Markov random fields (MRFs) are multivariate Gaussian distributions parameterised by a mean vector and a correlation matrix where there is a one-to-one relationship between the graphical representation of the model and the inverse correlation (precision) matrix [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF]. They are used in particular as the spatial random effects component within larger hierarchical models for the analysis of spatial data [START_REF] Banerjee | Hierarchical Modeling and Analysis for Spatial Data[END_REF][START_REF] Cressie | Statistics for Spatio-Temporal Data[END_REF].

We are interested in inference for the parameters of a Gaussian MRF with a spatial interpretation through the precision matrix. This inference is considered in a Bayesian framework where the model parameters are considered as random variables and the posterior density of the parameters is proportional to the likelihood of the model multiplied by a prior density of the parameters.

We first consider multivariate Gaussian distributions, their precision matrices and how this leads to a natural definition of Gaussian MRFs. Then we consider how Gaussian MRFs can be defined in terms of full conditional distributions which allows for the construction of the precision matrix of the whole model. We then outline finding the posterior densities for the model parameters which can be used for sampling and inference in a Bayesian framework.

Precision Matrix

In the following we consider vectors rather than sets of random variables but still use set notation in the subscripts, X V = [X 1 , X 2 , . . . , X n ] T when the meaning is unambiguous. For a matrix A we write the individual elements of the matrix as A ij and use standard vector and matrix notation.

X i |X V\i ∼ N j∈V\i ω ij X j , σ 2 i (3.1)
for all i ∈ V.

That is, each random variable has a conditional Gaussian distribution given all the other variables with mean, a weighted sum of the other variables and a given variance. We now consider finding the joint distribution of the CAR model.

Lemma 3.3 Brook's Lemma. Consider the collection of random variables X V with density

function p(x V ) > 0 for all x V ∈ Ω X V . For any x V , x V ∈ Ω X V we have p(x V ) p(x V ) = n i=1 p(x i |x 1 , . . . , x i-1 , x i+1 , . . . , x n ) p(x i |x 1 , . . . , x i-1 , x i+1 , . . . , x n ) (3.2) = n i=1 p(x i |x 1 , . . . , x i-1 , x i+1 , . . . , x n ) p(x i |x 1 , . . . , x i-1 , x i+1 , . . . , x n ) . (3.3)
Proof of Lemma 3.3. A proof is given by [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF].

Lemma 3.3 [START_REF] Brook | On the Distinction Between the Conditional Probability and the Joint Probability Approaches in the Specification of Nearest-Neighbour Systems[END_REF]) allows for the construction of the joint density of a CAR model from the full conditionals. However, there are conditions on the full conditionals in order that the density is proper. It is not the case that any arbitrary set of full conditionals will determine a proper distribution [START_REF] Gelfand | Proper Multivariate Conditional Autoregressive Models for Spatial Data Analysis[END_REF].

Theorem 3.4 Joint Density of a CAR Model. If the collection of random variables X V has full conditional distributions

X i |X V\i ∼ N φ j∈V\i e ij ∂ i X j , σ 2 ∂ i (3.4)
for all i ∈ V where e ij ≥ 0, e ij = e ji , e ii = 0, ∂ i = j∈V e ij and with parameters φ ∈ (-1, 1) and σ 2 > 0, the joint distribution of the model is

X V ∼ N (0, σ 2 (D(I -φW )) -1 ) (3.5)
where D = diag(∂ i ) and W is the symmetric matrix with entries e ij .

Proof of Theorem 3.4. First consider that for a zero-mean multivariate Gaussian distribution we have

log(p(x V )) ∝ - 1 2 x T V Qx V = - 1 2 i∈V j∈V Q ij x i x j = - 1 2 i∈V Q ii x 2 i + i∈V j∈V\i Q ij x i x j . (3.6)
Consider the full conditionals as given as (3.1). We have the density functions

p(x i |x V\i ) ∝ 1 σ i exp - 1 2σ 2 i (x i - j∈V\i ω ij x j ) 2
for all i ∈ V. Using (3.2) we can write

p(x V ) p(0) = n i=1 1 σi exp -1 2σ 2 i (x i - i-1 j=1 ω ij x j ) 2 1 σi exp -1 2σ 2 i ( i-1 j=1 ω ij x j ) 2 = n i=1 exp - 1 2σ 2 i x i - i-1 j=1 ω ij x j 2 + 1 2σ 2 i i-1 j=1 ω ij x j 2 = n i=1 exp - 1 2σ 2 i x 2 i -2 i-1 j=1 ω ij x i x j . Hence log p(x V ) p(0) = - 1 2 n i=1 1 σ 2 i x 2 i - n i=2 i-1 j=1 ω ij σ 2 i x i x j . (3.7)
Now we use (3.3) to obtain the corresponding expression

log p(x V ) p(0) = - 1 2 n i=1 1 σ 2 i x 2 i - n-1 i=1 n j=i+1 ω ij σ 2 i x i x j . (3.8)
Since (3.7) and (3.8) are equal,

log(p(x V )) ∝ - 1 2 n i=1 1 σ 2 i x 2 i + n i=1 j =i ω ij σ 2 i x i x j .
Comparing this to the form of the zero-mean multivariate Gaussian distribution (3.6), we have

Q ii = 1 σ 2 i and Q ij = ωij σ 2 i with ωij σ 2 i = ωji σ 2 j . Finally we write ω ij = φ eij ∂i and σ 2 i = σ 2 ∂i .
Hence the CAR model (3.5) is a Gaussian MRF with Q = 1 σ 2 D(I -φW ). There is a one-toone connection between the adjacency matrix W and the corresponding conditional independence graph since e ij = 0 if and only if Q ij = 0 and e ij > 0 if and only if (i, j) ∈ E. That is, the sum over j ∈ V\i in (3.4) may as well be over j ∈ ν(i) due to the zero pattern of W . Note that in general the adjacency matrix W is taken as a sparse matrix and hence Q is sparse. Recall that this is the point of constructing the precision matrix in this way and is necessary for efficient inference.

We additionally have interpretations for the parameters when considering X V as spatially embedded variables through the adjacency matrix W . If φ = 0 then the variables X V have a multivariate Gaussian distribution with a diagonal covariance matrix and hence the marginal distributions are independent. A value of φ > 0 results in spatial autocorrelation between neighbouring variables while a value of φ < 0 results in spatial anti-autocorrelation (easiest to see in the full conditional distributions (3.4)). Hence φ is known as the spatial autocorrelation parameter while σ 2 is a general variance term.

Posterior Densities

We consider inference for the parameters φ and σ 2 of the CAR model (3.5). Although CAR models are most often used as a spatial random effects component within larger hierarchical models in the analysis of spatial data [START_REF] Banerjee | Hierarchical Modeling and Analysis for Spatial Data[END_REF][START_REF] Cressie | Statistics for Spatio-Temporal Data[END_REF], inference for the parameters can also be considered when the models are used directly [START_REF] Bell | A Bayesian Analysis for Spatial Processes with Application to Disease Mapping[END_REF][START_REF] De Oliveira | Bayesian Analysis of Conditional Autoregressive Models[END_REF][START_REF] Ren | Objective Bayesian Analysis for CAR Models[END_REF]. This section is primarily based on the work of [START_REF] Bell | A Bayesian Analysis for Spatial Processes with Application to Disease Mapping[END_REF]. Theorem 3.5 Bayes' Theorem. For random variable Y with associated parameters θ and taking the parameters themselves as random variables we have

p(θ|y) = p(y|θ)p(θ) p(y) ∝ p(y|θ)p(θ).
Proof of Theorem 3.5. This is the definition of conditional probability density, p(x|z) = p(x,z) p(z) .

That is, the posterior density of the parameters is proportional to the likelihood of the model and a prior density of the parameters. This is the cornerstone of 'Bayesian inference' named after English mathematician Thomas Bayes . Rather than point estimates of model parameters, treating the parameters as random variables results in a posterior distribution that if not in an analytical form can usually be numerically sampled from.

Consider the CAR model with joint density given by (3.5) and parameters φ and σ 2 . We use the re-parameterisation σ 2 = τ -1 > 0 so that the model density can be written

p(x V |φ, τ ) ∝ |τ B(φ)| 1 2 exp - τ 2 x T V B(φ)x V
where B(φ) = D(I -φW ). We set priors p(φ) ∝ constant and p(τ ) ∝ τ -1 . We assume a priori that φ and τ are independent, that is p(φ, τ ) = p(φ)p(τ ) = τ -1 . Hence using Bayes' Theorem (Theorem 3.5) the joint posterior density is Now for fixed φ it can be seen from (3.9) that τ has a Gamma distribution. That is,

p(φ, τ |X V = x V ) ∝ p(x V |φ, τ )p(φ, τ ) ∝ |τ B(φ)| 1 2 exp - τ 2 x T V B(φ)x V τ -1 = |B(φ)| 1 2 τ n 2 -1 exp - τ 2 x T V B(φ)x V . ( 3
τ |{φ, X V = x V } ∼ Γ n 2 , 1 2 x T V B(φ)x V .
In order to find the posterior density of φ we integrate p(φ, τ |X V = x V ) with respect to τ using the identity (3.11) to obtain

p(φ|X V = x V ) = ∞ 0 p(φ, τ |X V = x V )dτ = |B(φ)| 1 2 ∞ 0 τ n 2 -1 exp - τ 2 x T V B(φ)x V dτ = |B(φ)| 1 2 1 2 x T V B(φ)x V -n 2 G n 2 ∝ |B(φ)| 1 2 (x T V B(φ)x V ) -n 2 .
(3.12)

Hence in order to sample from the joint posterior p(φ, τ

|X V = x V ) = p(τ |φ, X V = x V )p(φ|X V = x V )
we first sample from p(φ|X V = x V ) and then sample from p(τ |φ, X V = x V ) for fixed φ. ). By a regular grid we again mean that all vertices corresponding to the variables X V have the same degree (except the border vertices). This was an image analysis application concerning microscopy image data of cell culture models. In this case the vertices represent pixels in digital images and the modelling is carried out for each pixel and its four immediate neighbours, above, below, left and right in the digital image.

The graph for Publication 3 is also a hidden MRF but rather than a regular grid appears to be scale free. That is, in general terms, there are a few vertices with many edges and many vertices with a few edges. This isn't completely apparent in Figure 4.1 with only nine vertices constituting the Markov component of the model but we can still see the contrast to the regular grid case where every vertex has the same degree. Here the application was to do with genomic network data and the edges represent interactions between genes. While some genes may interact with most other genes, the majority of genes interact with only a few others.

The final graph corresponding to Publication 4 is also an MRF, however in this case the variables X V are not 'hidden'. In this application the vertices represent neurons and connections between vertices represent spatial interactions between neurons modelled by their relative spatial positions as obtained from image data. The neurons are not distributed within the image in the same way as the actual pixels, that is, they are not laid out in rows and columns. Hence each vertex may have a different number of edges but unlike the scale free case in Publication 3, the degree in this application is bounded in practice (see Publication 4). Note that final results for each application are summarised in Table 4.4 while the comparison and evaluation of our use of Markov models against the diverse literature associated to each application is covered within the publications themselves.

Publication 1 X a X b X c Z a Z b Z c Publication 3 X g X b X f X d X e X h X c X i X a Z g Z b Z f Z d Z e Z h Z c Z i Z a Publication 2 X g X h X i X d X e X f X a X b X c Z g Z h Z i Z d Z e Z f Z a Z b Z c Publication 4 X p X o X c X g X b X e X d X m X n X j X k X i X l X h X a X f
Our first application concerns time course gene expression data from grapevine plants undergoing berry development at multiple vineyards. The time series are modelled as HMMs where Z V are the observed gene expression measurements and X V are hidden underlying states (note that in Publication 1 the notation for hidden states is S and observed variables are X, W and C).

An alignment between expression profiles from different vineyards is achieved by parameterising a mapping to a common state sequence for each gene. Exact inference of the underlying 'hidden'

states is computationally feasible for HMMs owing to their tree structure while parameter inference is carried out with the EM algorithm [START_REF] Robinson | Alignment of Time Course Microarray Data with Hidden Markov Models[END_REF].

Both Publications 2 and 3 concern inference of the 'hidden' Markov component of hidden MRFs.

For Publication 2 this results in a segmentation of digital image data corresponding to different biological structures. The structures of interest are represented as contiguous objects within the digital images and hence using a MRF, where neighbouring pixels are induced to have the same segmentation label, is inherently appropriate for this application. Moreover, the smoothness of the segmentation output is particularly complimentary for the subsequent investigation of physical contact between the different cell types. The inference for Publication 2 is carried out using energy minimisation and requires the material from Chapter 2 concerning graph cuts and the α-expansion algorithm (Algorithm 2.2).

In Publication 3, genomic networks are similarly segmented by genes being labelled as a 'hit'

or not, that is, whether they are identified as valuable for further subsequent analysis. iments. Moreover we show that the alignment obtained under the model is robust against subsets of genes that do not align. We subsequently use a measure of alignment based on the inferred Viterbi paths to classify genes as likely to be developmentally driven or not and additionally validate the classification accuracy. This results in a set of ∼1200 genes with no current annotated function for which we have found new evidence that they are likely to be controlled in a developmental manner. tering allows for a segmentation with multiple labels to be achieved in an automatic and principled way. The segmentation output subsequently allows for the investigation of physical contact between tumour cells and CAFs on the surface of tumour organoids, which has important implications for tumour biology (Åkerfelt et al. 2015). We demonstrate how each step of our proposed method improves segmentation performance. Additionally we show our MRF based method is more generally applicable to other types of image data, not just fluorescence microscopy. proposed approach is inherently suitable for the complexity and resolution of the image data and allows for the analysis of upwards of thousands of neurons over hundreds of time points. We analyse data specifically developed for this study as well as from the literature. The quantification results both support and conform to previous analysis as well as demonstrate that our method has potential for screening applications.

Continued on the next page

Dans la publication 1 ,

 1 les modèles de Markov cachés (HMMs) sont utilisés pour atteindre un alignement et une classification des données de séries temporelles. Les publications 2 et 3 concernent l'inférence à la fois avec des MRFs cachés pour obtenir une segmentation des données d'images numériques et des données de réseau. Une analyse spatiale avec des MRFs de Gauss est présentée dans la publication 4. Nous montrons que notre utilisation particulière des modèles de Markov dans chaque publication nous a permis d'atteindre nos objectifs.

  Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Markov Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Hammersley-Clifford Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Inference for Hidden MRFs 2.1 Energy Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Graph Cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Alpha-Expansion Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Inference for Gaussian MRFs 3.1 Precision Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Conditional Autoregressive Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Posterior Densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 .

 1 Figure 1.1 shows conditional independence graphs corresponding to a number of different

Figure 1 . 1 :

 11 Figure 1.1: Conditional independence graphs corresponding to different Markov models.

Definition 2 . 3

 23 Directed Neighbours. In a directed graph -→ G , the directed neighbours of a vertex i are the members of the sets -

Algorithm 2 . 2

 22 Alpha-Expansion Algorithm. Given the collection of random variables X V with conditional independence graph G = (V, E), energy function (2.8), initial labelling x * V and specified label α ∈ L, 1. Construct the binary label problem for the dummy variables Y V as outlined above.

  .9) Definition 3.6 Gamma Distribution. The random variable Y has a Gamma distribution with parameters α > 0, β > 0 and support Ω Y = (0, ∞), written Y ∼ Γ(α, β), if and only if its

Figure 4 . 1 :

 41 Figure 4.1: Representative conditional independence graphs corresponding to each application.

  Each gene is assumed to have similar behaviour to the others it interacts with in the genomic experiments considered as represented by edges in the genomic network. We utilise an MRF to incorporate the additional network information in the identification of network based gene hits in the same way that adjacent pixels in a digital image are induced to have the same segmentation label. Hence Publication 3 also requires the material from Chapter 2 and we again use graph cuts and the α-expansion algorithm (Algorithm 2.2) for energy minimisation. Publication 4 concerns modelling the spatial component of neuronal synaptic activity in image data. In this case the variables X V are observed directly rather than 'hidden' and it is the spatial autocorrelation parameter within the Gaussian MRF (CAR) model that is of interest. Similar to the consideration of interacting genes, we model the mutual synaptic interaction of individual neurons through their neighbours as defined by relative spatial distance. Publication 4 requires the theory reviewed in Chapter 3 and sampling from the posterior distribution of the parameters of a CAR model in a Bayesian framework.

Publication 2

 2 We extend the established use of hidden Markov random fields (MRFs) for digital image segmentation of 3D fluorescence image data of cell culture models containing both multicellular tumour and cancer associated fibroblast (CAF) structures. Using both Gaussian mixture distributions and local entropy fil-

Table 4 . 4 -

 44 Continued from the previous page Publication 3 As an analogue to segmentation of image data, we use a hidden MRF based method to label 'hit' genes across genomic networks together with genomic screening data. The particular advance here concerns multiple 'hit' labels being easily achieved within the MRF framework. This allows for applications with multivariate measurements for each gene as well as further discrimination of the labelling. Such an approach had not been considered up until now and we demonstrate its advantages by finding additional pathway enrichment in previously analysed data concerning lymphoma. We show that our MRF based method is widely applicable to genomic screening data in general, including from RNA interference experiments. Additionally we compared to other competing methods using an independent simulation experiment and obtained the best performance results. Publication 4 We use a Gaussian MRF (CAR) model to quantify the synaptic activity phenotype of large populations of neurons. Previous analysis was limited to spike binning in an artificial square grid defined across the field of view whereas our model is based directly on the relative spatial positions of the individual neuronal somata. Spatial autocorrelation of the neuronal activity is then quantified through a Bayesian posterior density of the model parameters. Our

  Definition 1.1 Conditionally Independent Random Variables. † Random variables X and Y are conditionally independent given Z, written X ⊥ ⊥ Y | Z, if and only if

	p(x, y|z) = p(x|z)p(y|z)

  Equivalency of the Markov Properties. If the collection of random variables X V has density function p(x V ) > 0 for all x V ∈ Ω X V then the pairwise (1.1), global (1.2) and

	Theorem 1.3 local (1.3) Markov properties are equivalent.
	Proof of Theorem 1.3. The global Markov property clearly implies the local Markov property,
	which in turn implies the pairwise Markov property. From Lemma 1.2, all three properties are
	equivalent if the collection of random variables has a strictly positive density function.
	Hence a collection of random variables with strictly positive density function has a conditional
	independence graph for which Theorem 1.1 holds and hence is a Markov model. For all of the
	conditional independence graphs that we consider in the following we always assume that the
	corresponding collection of random variables has a strictly positive density function.
	Definition 1.7
	Definition 1.6 is given in terms of the local Markov property (1.3) which is the classic definition
	of a Markov chain. That is, each random variable is conditionally independent of all other random
	variables in the model given its 'neighbours' as represented in the conditional independence graph.

  Definition 2.7 Maximum Flow and Minimum Cut. A cut [S] is a minimum cut if there exists no other cut [S * ] such that val[S * ] < val[S]. A flow f is a maximum flow if there exists no other flow f * such that val(f ) < val(f * ). Theorem 2.2 Maximum Flow/Minimum Cut. For a network N , the value of the maximum flow is equal to the value of the minimum cut.

  Definition 3.5 Conditional Autoregressive Model. The collection of random variables X V is a conditional autoregressive (CAR) model with parameters ω ij ∈ R and σ 2 i > 0 if and only if the collection of full conditionals distributions are

Definition 3.1 Symmetric Matrix. A square matrix A is symmetric if and only if A T = A.

Table 4 .

 4 2:The graphical model structure corresponding to each application. chain for Publication 1 we mean a classic Markov chain where all vertices corresponding to the variables X V have the same degree (except the beginning and end vertices). This application concerned time course gene expression data and we model each instance of the gene expression in time using a hidden Markov model (HMM). In Publication 2 we are now considering a non-tree Markov random field (MRF

		Non-tree Variables X V are 'hidden' Markov Model	Description
	Publication 1		X	HMM	Regular chain
	Publication 2	X	X	Hidden MRF	Regular grid
	Publication 3	X	X	Hidden MRF	Scale free
	Publication 4	X		Gaussian MRF Triangulation
	By a regular

Table 4 .

 4 3: The inference carried out in each application. We are now able to consider how each of the Markov models are utilised, what kind of inference is carried out and for what aims based on the type of Markov model and what the vertices and edges represent in each application. As noted, the most important information in Table4.2 is given in the first two columns corresponding to the particular type of Markov model. This determines the inference that can be achieved in each application, which is set out in Table4.3 and discussed below.

		Infer 'hidden' X V Infer parameters	Framework	Reference
	Publication 1	X	X	Maximum	(Robinson 2012)
				likelihood
	Publication 2	X	X	Maximum	Chapter 2
			X	a posteriori
	Publication 3	X		Maximum	Chapter 2
				a posteriori
	Publication 4	X	X	Bayesian	Chapter 3
		X

Table 4 .

 4 [START_REF] Robinson | Spatial Quantification of the Synaptic Activity Phenotype Across Large Populations of Neurons with Markov Random Fields[END_REF] summarises the results for each application. Hence we have shown that Markov models can be used to generate new knowledge in diverse applications relating to multiple domains of biological research. The Markov structure allows for a graphical and intuitive model of conditional independencies within a standard and general framework, making such models very popular in many different areas. Furthermore, as we have shown, there are many associated computational schemes and algorithms available in order to efficiently utilise Markov models in practice.

Table 4 .

 4 4: The results for each application. Publication 1 We propose a novel alignment method for time series gene expression data based on hidden Markov models (HMMs). The modelling extensions we employ account for specific features of the data, in particular the sparsity of the time series commonly encountered in such gene expression time course exper-

Acknowledgements

am truly grateful to Anni for all her love and support which has meant so much to me.

Funding for the first three years was jointly provided by VTT, Finland and CEA, France while finding for the final fourth year was provided by the MATTI doctoral school at the University of Turku and the Finnish Cultural Foundation, Varsinais-Suomi Regional Fund. Funding to participate in a number of international conferences as well as travel between Finland and France was provided by both labs at VTT and CEA, the Turku University Foundation, the MATTI doctoral

Proof of Lemma 3.1. Symmetry is clear and for any vector u = 0 take v = A -1 u. Since A is invertible, v = 0 and hence v T Av > 0 implies u T A -1 u > 0.

Definition 3.3 Multivariate Gaussian Distribution. The collection of random variables X V = [X 1 , X 2 , . . . , X n ] T has an n-variate Gaussian distribution with mean vector µ ∈ R n , symmetric positive definite covariance matrix Σ ∈ R n×n and support Ω X V = R n , written X V ∼ N (µ, Σ), if and only if its density function has the form

Gaussian distributions are named after German mathematician Carl Friedrich Gauss (1777-1855). Consider the matrix Q = Σ -1 known as the precision matrix and note that since Σ is symmetric positive definite, Q is symmetric positive definite (Lemma 3.1). There is a direct connection between the Markov properties of the model and the precision matrix Q. The content of this section is based on the presentation of [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF].

Theorem 3.2 Conditional Independence and the Precision Matrix. Let the collection of random variables X V have a multivariate Gaussian distribution with precision matrix Q = Σ -1 .

Then

for all i = j ∈ V.

Proof of Theorem 3.2. A proof is given by [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF].

That is, for any multivariate Gaussian model, we have a connection between the form of the precision matrix and the conditional independence properties of the model. The covariance matrix Σ gives information about the marginal dependence structure of the model whereas the precision matrix Q gives information about the conditional independence structure of the model. This leads to a natural definition of a Gaussian MRF through the precision matrix.

Definition 3.4 Gaussian MRF. The collection of random variables X V is a Gaussian MRF with respect to a graph G = (V, E) with mean vector µ ∈ R n , symmetric positive definite precision matrix Q ∈ R n×n and support Ω X V = R n if and only if its density function has the form

and

It is clear that the three Markov properties (pairwise (1.1), global (1.2) and local (1.3)) are all equivalent for Gaussian MRFs since p(x V ) > 0. Hence a Gaussian MRF is defined through the precision matrix Q which is one-to-one with a graphical representation of the model. The non-zero

no edge between vertices i and j in the conditional independence graph and hence X i and X j are conditionally independent given X V\{i,j} . So we can simply read from the elements of Q whether two variables in the model are conditionally independent. Now Q is generally taken as a sparse matrix (most of its elements are zero) while Σ may be very dense (most elements are non-zero). The Markov structure of a Gaussian MRF allows for efficient inference due to a sparse precision matrix, which corresponds to a sparse conditional independence graph (Chapter 1). Although it is possible to define reasonable dependence structures through a covariance function, there are no guarantees on the sparsity of the inverse covariance matrix and it is non-trivial to construct reasonable dependence structures ensuring sparsity [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF].

So rather than considering covariance matrices, we consider conditional independencies and the construction of reasonable and useful precision matrices.

Conditional Autoregressive Models

To utilise a Gaussian MRF in practice, instead of specifying the full covariance matrix for every pair of variables, we consider the 'local' behaviour for each variable and the full conditional distributions, as introduced by [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems (with Discussion)[END_REF]. That is, rather than defining a joint density p(x V ), we consider the much simpler task of defining the conditional densities p(x i |x V\i ) for all i ∈ V.

We consider some final details about sampling using the posterior density function (3.12).

Firstly, in applications with large n, the calculation of the determinant of B(φ) can easily become computationally infeasible while the exponentn 2 will result in computational zeros. For the later problem we consider log space (scaling and then converting back as necessary) and write

Hence the exponent is no longer a problem but a potential problem with the calculation of the determinant remains.

Definition 3.7 Cholesky Decomposition. Every symmetric positive definite matrix A can be written as

where L is a unique lower triangular matrix with positive diagonal entries.

We use this decomposition, named after French mathematician André-Louis Cholesky , to rewrite ) is greatly simplified. Although finding the Cholesky decomposition may be difficult in general, since Q is sparse it is possible to use numerical methods for sparse matrices to efficiently find the Cholesky decomposition and hence use CAR models in practice [START_REF] Rue | Gaussian Markov Random Fields: Theory and Applications[END_REF].

Chapter 4

Presentation of Publications

We make use of Markov models in four bioinformatics applications each of which has a corresponding publication presented below. Firstly, Table 4. 4.2. Note that in this discussion we focus on the Markov components X V within each model. The most important differences between the models as seen in Figure 4.1 are whether the models are trees or not and whether the variables X V are 'hidden'

or not (compare to Figure 1.1). Recall that these difference define the particular type of Markov model and determine the possibilities for inference. We briefly consider general descriptions of the models and how they relate to each application.

Following my Master of Philosophy thesis [START_REF] Robinson | Alignment of Time Course Microarray Data with Hidden Markov Models[END_REF], the underlying research framework of the use of Markov models was driven by me while biologically relevant applications were developed within multiple research groups. Table 4.5 lists the published author contributions.

For each application I was primarily responsible for:

• Methodological development, implementation and analysis of results

• Literature review, evaluation and presentation of results

• Production of figures, movies and code for publication

• Writing the manuscript and organising the contributions of co-authors

• Submission of the manuscript and responding to reviewers Table 4.5: The author contributions for each application.


Laurent Cassayre

Torstein A Utigard

Sylvie Bouvet

Visualizing Gas Evolution on Graphite and Oxygen-Evolving Anodes

Recent progress in material science might soon allow the replacement of the consumable carbon anode by an inert material. This is likely to induce changes in the overall process, and particularly in the gas evolution. Video recordings of oxygen-evolving anodes (SnO 2 , Cu, Cu-Ni) and carbon anodes were performed in laboratory electrolysis cells, using direct observation from above, a see-through cell, and radiography techniques. The gas behavior was very different between the two kinds of anodes, and probably linked to the wettability of the material by the electrolyte.

INTRODUCTION

Replacement of carbon anodes by non-consumable (or inert) materials has been a subject of research since the invention of the Hall-Héroult process more than a hundred years ago. With the proper material, inert anodes would offer potential savings in operating costs by eliminating carbon anode replacement and by providing tighter control of the anode-cathode distance. As presented by the following equations, the use of an oxygen-evolving anode (Equation 1) could also help suppress the emissions of greenhouse gases, which are a growing concern with the current process (Equation 2).

Al 2 O 3 (dissolved) → 2 Al (l) + 3/2 O 2 (g) (1) Al 2 O 3 (dissolved) + 3/2 C (s) → 2 Al (l) + 3/2 CO 2 (g) (2) While a vast array of materials has already been patented for inert-anode application, none is yet in commercial use. This is mainly because the severe operating environment and the need for good electronic conductivity limits the number of suitable materials. Most recent research has concentrated on metals, oxide ceramics, and cermets. [START_REF] Pawlek | Inert Anodes for the Primary Aluminium Industry: An Update[END_REF] In the present study, copper, coppernickel alloy (wt.% 75/25), and tinoxide materials were tested as oxygenevolving anodes, even if none of them is stable enough to totally resist corrosion.

Modifying the anodic material is likely to induce many changes in the overall process; cell voltage, bath resistivity, and gas evolution are, indeed, deeply linked with anode properties. The aim of the current work was to compare the gas behavior between oxygen-evolving and graphite anodes, and to characterize the bubble size and the bubble-layer depth in the interelectrode area. Video recordings of electrolysis at various current densities (CD) were performed in laboratory cells, using direct observation from above and see-through-cell 2-4 and radiography [START_REF] Utigard | Direct Observation of the Anode Effect by Radiography[END_REF][START_REF] Utigard | Visualization of the Hall-Héroult Process[END_REF] techniques. The influence of the alumina content in the electrolyte was studied for graphite and tin oxide. A comparison of the wettability of the materials by the electrolyte was also carried out to try to explain the very different gas behavior observed.

EXPERIMENTAL SETTINGS

Three different methods were used to get video recordings of the gas evolution at the anode during alumina electrolysis. Schematics of the experimental techniques are presented in Figure 1.

Direct observations from above the bath were carried out in a small furnace open to the air with a microscope positioned above the cell. The cathode was a graphite cylinder (∅ = 5 mm), while the anode was a graphite or tin-oxide cylinder (∅ = 5 mm).

Electrolyte

Electrolysis was also carried out in transparent quartz tubes (∅ in = 54 mm) placed in an electric furnace, allowing lateral observations through a hole in the front panel. In this see-through cell, the cathode was a graphite cylinder (∅ = 43 mm, h = 15 mm) located at the bottom of the tube, and the anode was either a graphite or tin-oxide cylinder (∅ = 9 mm, h = 50 mm), vertically positioned 2 cm above the cathode. A new quartz tube was used for each experiment.

Radiography tests consisted of performing electrolysis in a graphite crucible (∅ in = 45 mm) placed in a vertical quartz tube and heated in a di-moly silicide resistor furnace. Images of the anode immersed in the electrolyte were obtained from an x-ray beam passing through whole system. An electrically insulating alumina cylinder protected the inside walls of the crucible from the electrolyte so that current distribution was kept vertical through the bath. A boron-nitride cylinder was used for the same purpose with low alumina-content electrolytes. Anodes (∅ = 18 mm, h = 50 mm) were made of ATJ graphite, tin oxide, pure copper, or copper-nickel alloy (wt.% 75/25). The atmosphere in the cell was purged with nitrogen before heating, and a small flow was kept during the runs to prevent carbon oxidation.

Two different electrolyte compositions were used, based on a mix of cryolite with AlF 3 (11 wt.%), Al 2 O 3 , and CaF 2 (5 wt.%). Weight fraction of alumina was either 2.5 wt.% (electrolyte A) or 9 wt.% (electrolyte B), electrolyte B being saturated in alumina. [START_REF] Solheim | Liquidus Temperature and Alumina Solubility in the System Na 3 AlF 6 -AlF 3 -LiF-CaF 2 -MgF 2[END_REF] The mix of powders was pre-melted at 1,000°C in graphite crucibles, then cast, crushed, and stored in a dessicator.

BUBBLING AT THE GRAPHITE ANODE

Bubble formation on the side of a graphite anode was observed from above at small CDs (0.05 to 0.2 A.cm -2 ). Bubble birth occurs on specific nucleation spots, at a frequency growing with CD. Increasing CD leads also to a large number of nucleation sites. Roundshaped bubbles grow and coalesce until reaching 2-3 mm in diameter, then detach and escape vertically. As shown in Figure 2, the surface coverage of the anode is very high. The bubbles are stuck on the anode and, without much agitation in the bath, detach with difficulty.

Tests in see-through cells were limited in time (no more than 7 minutes) and CD (up to 0.5 A⋅cm -2 ), because streamers rapidly arose in the cell and lead to an opaque electrolyte. This fog was attributed to dissolved aluminum produced at the cathode, [START_REF] Haupin | See-through Hall-Héroult Cell[END_REF][START_REF] Zhuxian | Formation of Metal Fog during Molten Salt Electrolysis Observed in See-through Cell[END_REF] and did not fade as reported by some authors. [START_REF] Zhuxian | [END_REF] Observations at industrial CDs (around 0.7 A⋅cm -2 ) could not be performed, diminishing the interest of the method.

Video recordings of gas evolution during both see-through-cell and radiographic experiments confirmed previous observations described in the literature: [START_REF] Xue | Bubble Behaviour-Cell Voltage Oscillation during Aluminium Electrolysis and the Effects of Sound and Ultrasound[END_REF][START_REF] Utigard | Anode Gas Behaviour during Electrolysis[END_REF] a periodic gas mechanism release occurs at the bottom of the anode, leading to important variations with time in the amount of gas trapped in the interelectrode area. At small CDs (up to around 0.5 A⋅cm -2 ), several bubbles form on the bottom surface of the anode, then grow and coalesce until a single large bubble covers the whole area and is finally released, leaving a clean surface free of gas. Figure 3 shows the bubble trapped under the anode, just before release, during see-through cell electrolysis at 0.2 A⋅cm -2 .

At higher CDs (studied only by x-ray experiments), coalescence does not happen as much as at low CDs, and bubbles escape before covering the whole anode, as illustrated in Figure 4.

The average bubble diameter before release, calculated as the mean diameter of ten bubbles for each CD, is presented in Figure 5. The diameter decreases from around 18 mm at 0.2 A⋅cm -2 to 7 mm at 1.6 A⋅cm -2 . The alumina content in the electrolyte does not have a major influence on the bubble size.

The release frequency of the CO 2 bubbles is also strongly related to the current density, as shown in Figure 6. These values agree with those obtained on 10 cm diameter anodes by monitoring the fluctuations of the bath level, 10 but significantly higher than those calculated by cell voltage analysis 3 using 1 cm diameter anodes.

BUBBLING AT OXYGEN-EVOLVING ANODES

Gas evolution is nearly identical for tin-oxide, copper, and copper-nickel anodes, but was found to be very different as compared to what was observed on the graphite anode.

Observations of a tin-oxide anode from above showed tiny bubbles escaping continuously without much coalescence, forming a foam around the electrode, as illustrated in Figure 7. The bubble diameter, measured under the microscope (40 x), was estimated to be around 0.1 mm, which is ten to 30 times smaller than on graphite anodes.

Some tests were run in see-through cells, but observations were difficult because the froth of bubbles hid the electrode even at low CDs (0.1 A⋅cm -2 ). This resulted in poor contrast between the electrolyte and the gas. Furthermore, bath mist occurred very quickly, the same way as described when using graphite anodes.

Radiographic observation proved to be more adapted to the study of the gas evolution on oxygen-evolving anodes. The anode appears in black and the bath in dark grey, while gas takes the form of a light grey halo around and under the electrode (Figure 8). As almost no coalescence occurs between the tiny bubbles, gas release looks continuous and the gas foam under the electrode keeps roughly a constant shape at fixed CD.

The depth of the gas layer under the anode is plotted as a function of CD and material in Figure 9. The evolution is almost identical for the three oxygenevolving electrodes: the bubble layer thickness grows regularly with CD until reaching 2-2.5 mm at 1 A⋅cm -2 . Increasing CD over 1 A⋅cm -2 doesn't increase the penetration of the gas under the anode, maybe because the volume fraction of the gas in the bubble layer increases. However, the contrast of the x-ray images was not sufficient to observe noticeable changes in the brightness of the gas froth, as it could be expected with a higher gas density. Another possibility is that the release velocity of the bubbles increases, leading to a shorter residence time of the bubbles under the anode.

When testing electrolyte A (2.5% Al 2 O 3 ) on a tin-oxide anode, no difference was seen in the gas behavior as compared to the electrolyte B (9% Al 2 O 3 ).

Figure 9 shows that the bubble layer under the graphite anode is decreasing with CD, from 5 mm to 4.2 mm, confirming both the values and the trend reported by Aaberg et al. [START_REF] Aaberg | The Gas under Anodes in Aluminium Smelting Cells Part II: Gas Volume and Bubble Layer Characteristics[END_REF] At 1 A⋅cm -2 , gas penetration under the anode is roughly twice as deep with graphite than with oxygen-evolving anodes, which is important with regards to the setting of the anode-cathode distance.

WETTING OF THE ANODES

The angle of contact, θ, is a measure of the degree of wetting of a solid surface by a liquid phase. In the case of regular aluminum electrolysis, θ is dependant on the composition of the melt, the physical properties of the carbon material, the temperature, the composition of the gaseous phase, and the current density. [START_REF] Grjotheim | Aluminium Electrolysis: Fundamentals of the Hall-Héroult Process[END_REF] An interesting point is that, as suggested by Vogt, 12 a decrease of this contact angle might be correlated to an increased bubble size. It was then decided to use a hot-stage microscope to study the wetting of graphite, tin-oxide, copper, and coppernickel alloy by the electrolyte, in order to explain the difference in bubble size observed during electrolysis.

Thirty-five milligrams of electrolyte B (9% Al 2 O 3 ) were placed on a horizontal solid sample (∅ = 10 mm) and quickly heated to around 1,000°C. A video camera mounted on a microscope recorded the melting of the bath on the sample. Experiments were usually carried out in the air, except for some tests on graphite samples performed in CO 2 .

The contact angle of molten electrolyte on a graphite sample was measured around 120-130°, which is consistent with literature values. [START_REF] Grjotheim | Aluminium Electrolysis: Fundamentals of the Hall-Héroult Process[END_REF] The influence of the gaseous phase (air or CO 2 ) on the wettability could not be determined; it was not significant compared to the measurement uncertainties of the contact angle (±5°). Graphite is thus poorly wetted by the electrolyte, as illustrated in Figure 10.

Wetting on tin oxide was found to be very different. As shown in Figure 11, the electrolyte started to spread as soon as it began to melt, and then rapidly covered the top surface and the sides of the sample. Contact angle is therefore close to zero. Copper and copper-nickel samples were also very well wetted by the electrolyte.

According to the preceding results, the size of the bubbles released during alumina electrolysis could be related to the wettability of the anode by the electrolyte. The bubbles formed on tin oxide, copper, and copper-nickel, all very well wetted by the electrolyte, will detach easily from their nucleation site, while bubbles formed on carbon are stuck and grow before detachment.

CONCLUSIONS

The characterization of the bubble layer under the anode, summarized in Table I, showed many differences between the two types of anodes, which could have important consequences in the overall process:

• Given that the bubble layer is thin and keeps a constant shape on oxygen-evolving anodes, a smaller anode-cathode distance than currently used on carbon anodes could be set and help save energy by reducing the ohmic losses through the bath. • Bath stirring, which is necessary to keep a constant alumina concentration close to the electrodes, might be significantly reduced. • The back reaction (reoxidation of the aluminum produced at the cathode and dissolved in the bath) might be favored by the increase of the gas/electrolyte area due to the small size of the bubbles. The study of three different oxygenevolving anodes showed that a good wetting by the electrolyte leads to the formation of tiny bubbles. Contact angle measurements of the system electrolyte-anode-oxygen might then allow predicting the gas evolution on other candidate as inert anodes materials without having recourse to heavy equipment like the x-ray furnace.

Figure 1 .

 1 Figure 1. Schematics of the experimental conditions: (a) Observation from above; (b) see-through cell; and (c) x-ray furnace.

Microscope

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. An observation from above-side bubbles on graphite anode. [∅ anode = 5 mm, CD = 0.2 A⋅cm -2 , T = 960°C, electrolyte B (9% Al 2 O 3 )]

Figure 6 .Figure 5 .Figure 7 .

 657 Figure 6. The evolution of bubble release frequency with CD under a graphite anode. [∅ anode = 28 mm, T = 960°C, electrolyte A (2.5% Al 2 O 3 ), ACD = 2 cm]

Figure 10 .Figure 8 .Figure 9 .∅

 1089 Figure 10. Molten electrolyte on graphite at 1,000°C.

ACKNOWLEDGEMENTS

The authors are grateful to Pechiney for the financial support of this work.