Summary

• Heart diseases

• Inflammatory diseases

• Brain tumours

  ---

• Image registration tool

•  NEWGenerative model of heart growth using Exoshape NEW

Contributors

Heart diseases

Integrated Models of Paediatric Heart Diseases

I- Knowledge-Based Models of Paediatric Heart Diseases

Measurement of the right ventricle volume from
ultrasonography images. (Image from Gaslini)

Abstract. During the diagnosticprocess, cardiologists and clinicians perform several measurements toidentify the pathologies and to assess their severity. One of the mainpurposes of Health-e-Child is to gather the knowledge and the experienceof the cardiologists involved in the project and to designdecision-support models that can help the diagnostic of paediatric heartpathologies. In this section we present the underlying principles of apreliminary knowledge-based model that is currently in development. Itmainly relies on a set of representative measurements and on a modelthat formalises the clinical methods used during patients evaluation toautomatically identify and classify heart pathologies.

1. Knowledge-based models of paediatric heart diseases

II- Disease-Based Models of Cardiac Electromechanics

AHA segmental division of a generic anatomy

Abstract. In the last decades thecomputational modelling of the human body, and in particular of theheart, has attracted increasing amounts of interest (Ayache, 2004) as it hasbenefited from progresses in biology, physics and computer science. Theunderstanding of the heart functions, from the heart anatomy to theelectrical activity to the biomechanical phenomena, is crucial for thedevelopment of advanced treatments of heart diseases. We describe herean in-silico 3-D model for the simulation of the heart functionin healthy and ill children. Based on the approach proposed by(Sermesant et al., 2006), this model includes knowledge coming fromvarious disciplines (anatomy, electrophysiology and biomechanics) in aframework that allows direct evaluation with respect to actual clinicalmeasurements.

The model is first adjusted to thephysiology of the normal heart in order to calibrate theparameters. Then, preliminary simulations of the paediatric cardiacdiseases considered in Health-e-Child are performed by modifying theappropriate macroscopic parameters. Although basic, these simulationsprovide promising results, the clinical parameters vary accordinglywith respect to the simulated pathology.

1. Electromechanical model of the heart
2. Simulation of the normal heart
3. Simulation of right-ventricle overload
4. Simulation of dilated cardiomyopathy
5. Simulation of hypertrophic cardiomyopathy
6. Conclusion

III- Patient-Specific Anatomical Models of the Heart for Decision Support System

Segmentation of the heart from 3-D MRI.

Abstract. In the previous sectionstwo different methods have been described for the assessment and thesimulation of paediatric heart diseases. On one hand, there are themodels that rely exclusively on parameters representing the patient'sphysiology; and on the other hand, there are the models where onlygeneral considerations are used. Albeit efficient, these two methods donot take advantage of all the available information. The understandingof paediatric cardiac pathologies would greatly benefit from thecombination of the two approaches.

An effort must then be done to adapt to thepatient physiology the general framework described inthe second section of thiswebsite. A common approach consists in using 3-D imaging andsegmentation techniques to directly create a 3-D representation of thepatient's anatomy. However, these methods are often limited by thequality of the underlying image. Another strategy would be to build ageometrical model, as in the disease-based models, but by adapting thesynthetic geometry to clinical parameters measured during clinicalevaluation. Such an approach would allow us to get a good approximationof the patient's heart anatomy without the requirement of anyhigh-resolution 3-D image.

1. Current method: From MRI segmentation to 3-D anatomical model of the heart
2. Perspective: From 1-D clinical parameters to 3-D anatomical model of the heart

IV- Personalised Simulation of Myocardium Electromechanics and Pulmonary Valve Replacement Surgery in Repaired Tetralogy of Fallot: a Case Study

Abstract.So far, we have presented on the one hand an electromechanical model of the heart, which is able to perform generic simulation of congenital heart diseases, and on the other hand we described methods to segment the heart from medical images. We investigate in this section how we can use both tools to perform personalised simulations of the heart, in a Health-e-Child patient, and simulate therapies.

We focus on repaired Tetralogy of Fallot. For these patients, pulmonary valve replacement (PVR) is an effective treatment. The therapy consists in replacing the pulmonary valves, surgically or through minimally invasive approaches. The surgeon can also directly remodel the right ventricle to remove any impaired tissue. However, the postoperative outcomes of these techniques mostly depend on patient patho-physiology. We propose in this section an integrated workflow to simulate and predict the effects of PVR upon the cardiac function of a patient. We first personalise an active, anisotropic, visco-elastic electromechanical model of the heart. The two different interventions are then simulated. PPVR is achieved by disabling the regurgitations. Volume reduction surgery is simulated using a soft-tissue intervention platform. Tested on a young patient, the framework showed promising results. Preoperative simulations were successfully personalised (reference RV ejection fraction (EF): 41%, simulated EF: 40%). Postoperative simulations showed that, on this specific patient, the best improvement in cardiac function was achieved by PVR with surgery (EF after PPVR: 40%, after PVR with surgery: 51%).

1. Introduction: Medical Background and Challenges
2. Creating a Personalised Anatomical Model of the Heart
3. Personalising an Electromechanical Model to Simulate Patient Cardiac Function
4. Simulating Pulmonary Valve Replacement Therapies
5. Conclusions

NEW Automatic Estimation of 3D Myocardium Strain from Anatomical Cardiac Images NEW

Automatic estimation of myocardium deformation
(Data courtesy King's College London, St Thomas Hospital,
Division of Imaging Sciences)

Abstract.The heart is a moving organ and evaluating its intrinsic deformationis crucial for a comprehensive evaluation of its condition. Imagingmodalities have been developed to measure the cardiac deformation butwhen these technologies are not available, little can be provided tothe clinicians to assist his diagnostic. Unfortunately, this is thetypical situation in congenital heart diseases like tetralogy ofFallot, where only echocardiography or standard cMRI are available sofar. In Health-e-Child we tackle the objective of estimating thecardiac deformation from dynamic anatomical images like standard cineMRI. To that end, we constrain a demons-based non-linear imageregistration algorithm to provide incompressible deformations, a majorfeature of the myocardium, to cope with the scarce textures and thelow image quality. This leads us to a well-posed framework with newinsights on demons regularisation. This section describes the overallmethodology and an application on a real case. Reported results arepublished in (Mansi et al., 2010)

1. Introduction and Method: iLogDemons, a Consistent Incompressible Image Registration Algorithm
2. Experiments on Synthetic Data
3. Experiment on a Real Case

References

• Ayache, N., Ed., 2004. Computational Models for the Human Body, ser. Handbook of Numerical Analysis, P. Ciarlet, Ed. Amsterdam, The Netherlands: Elsevier.
• Sermesant, M., Delingette, H., Ayache, N., 2006. An Electromechanical Model of the Heart for Image Analysis and Simulation. IEEE Trans. Med. Imag., vol. 25, no. 5, pp. 612-625.
• T. Mansi, X. Pennec, M. Sermesant, H. Delingette, and N. Ayache. LogDemons Revisited: Consistent Regularisation and Incompressibility Constraint for Soft Tissue Tracking in Medical Images. In Medical Image Computing and Computer Assisted Intervention (MICCAI), Lecture Notes in Computer Science. Springer, 2010. In press.

Knowledge-Based Models of Paediatric Heart Diseases

1. Introduction

During clinical evaluation, cardiologists and clinicians measure several key parameters to identify the pathologies and to assess their severity. Ventricle volumes, contractility and blood flows for example are essential features of the cardiac function, and abnormalities in their values can be caused by heart diseases.

One of the main purposes of Health-e-Child is to gather the knowledge and the experience of the cardiologists involved in the project and to design a decision-support model that can help the community to detect and assess paediatric heart pathologies. This section is the knowledge-based model that results from the analysis of the methods used by cardiologists during the diagnostic process. It is based on the interpretation and the classification of some cardiac parameters by means of left-/right-ventricular comparisons and basic group analyses.

2. Estimation of cardiac parameters

Ultrasonography and MR imaging allow the clinicians to measure geometrical and functional parameters of the heart. From these data, they can estimate for example the volume and mass of either the whole heart or segments of interest, the blood flows, etc. More subtle measurements such as the global or segmental contractility of the myocardium can also be derived from these images, provided that 2D+t or 3D+t acquisitions are available.

In Health-e-Child, the volume of the cavities, the muscle contractility and the blood flows across the cardiac valves are estimated from 2-D/3-D ultrasonography and cardiac MR imaging.

To acquire the total volume of the right ventricle, advanced 3-D echocardiography sequences tailored to the clinical evaluation of the right ventricle (RV) in children are used (Gaslini). On the other hand, because such modality is not available for the evaluation of the left ventricle, the total volume of this cavity is estimated by using MR imaging only. Of course, when it is possible, RV volume is also measured through MRI, improving in this way previous measurements performed with ultrasonography. In both cases, the recommendation of the American Society of Echocardiography are complied.

Segmental volumes are also estimated to detect local abnormalities such as anomalous wall motions. Right ventricle for example is commonly subdivided into three segments: inlet, outlet and apical-trabecular (see next video).

Measurement of the volume of the right ventricle from ultrasonography images. (Image from Gaslini)

Get the Flash Player to see this player.

Volumetric representation of the right ventricle. Yellow: outlet, Green: inlet, Red: apical. (Video from Gaslini)

The estimation of the volumes of the heart cavities from MR imaging is performed by using the modified Simpson's algorithm. Contrary to 2-D ultrasonography evaluations and to other algorithms, the Simpson's method does not rely on any assumptions about the heart geometry. The ventricles are first divided into 12-15 slices of approximately 8mm. The areas of the cavities are then measured at each slice and, by multiplying them by the slice thickness, the intermediate "truncated" volumes are obtained. The total volume is finally calculated by summing all the intermediate results. Owing to this algorithm and to a better image resolution, MRI clinical evaluations are often considered more reliable and more accurate.

Estimation of the volume of the ventricles on MRI.

The second set of parameters consists of blood flow measurements. During the cardiac cycle, the blood can go from the atrium to the ventricle during diastole and from the ventricle to the artery during systole. The estimation of the amount of blood that crosses the valves can thus provide useful information about the heart function. Anomalous values may point out pathological features of the myocardium, ventricular overload or valvar regurgitations. For instance, tricuspid flows that are too high with respect to the mitral flows may suggest the presence of shunt between the atria.

Flow measurements are thus crucial for diagnosis. To estimate them, ultrasonography imaging is usually used and in particular the Doppler modality. This technology relies on the Doppler effect to estimate the velocity field of the observed structures.

Concretely, blood flows across the valves are assessed as follows. The specialist first estimates the area of the valve by using ultrasonography. Then, he gets the flow velocity at the level of the valve of interest, during several cardiac cycles, by means of continuous-wave Doppler acquisitions. Finally, he obtains the mean flow that passes across this valve by using the following formula:

The following video illustrates the Doppler acquisition process. At each time t, the expert can visualise the blood velocity and the 2D plane under study. The flow that crosses the white dot line in the US image corresponds to the area of the region bounded by the velocity curve and the x-axis.

The last parameters that are acquired for the knowledge-based models are contractility measurements. They quantify the efficiency of the cardiac muscle. As for the previous data, those parameters are mainly obtained by using ultrasonography images and, if available, from cine MRI (3D+t).

A global contractility index is measured by considering the entire ventricle, while segmental contractilities can be acquired for each segment of the heart. The ejection fraction of each ventricle (the amount of blood in percentage ejected during a cardiac cycle) is also estimated to evaluate the efficiency of the organ.

To conclude, all the clinical measurements that can be used in the knowledge-based model are summarised in the following table.

3. Pre-processing of the acquired data

Once acquired, the raw clinical parameters must be pre-processed before doing any group analyses or comparisons. In this way, the values are first normalised, followed by the calculation of standard scores to point out significant differences between the patient and a population of reference.

The first step of the pre-processing stage aims at normalising the acquired values in order to remove the possible effects of the height and weight of the patient. Indeed, height and weight can have a great influence upon the cardiac function and the acquired data must therefore be normalised to cancel their impact.

In Health-e-Child, the clinical measurements are normalised by dividing them by the body surface index computed by using Dubois' equation:

After normalisation, standard scores (also known as Z-score) are computed. The standard score is a dimensionless quantity that reveals how many units of the standard deviation the patient's parameters differ from the normal values (estimated among a population of normal subjects). It is calculated by using the formula:

Values higher or lower than two standard deviations (z-score higher or lower than 2) are considered as significant and may point out abnormal cardiac function.

4. Simplified Description of the Clinical Workflow

The imaging data presented in the previous section are integrated and used as input to a decision process. An exemplar system could be: the parameters are compared with normal values and displayed to the cardiologist in a user-friendly way; ventricular segments with abnormally high or low contractility are highlighted and measurements greater or lower than 2 standard deviations are pointed out; finally, Q_p/Q_s ratio is displayed.

In addition to the direct visualisation of the acquired parameters, right-ventricle measurements are compared with those of the left ventricle. This comparison relies on the fact that the circulatory system is a closed system and therefore the ratio of any parameter of the right ventricle over the equivalent parameter of the left ventricle must be equal to 1±5%. An anomaly in these ratios may be an indicator of pathological features.

First step: flow comparisons

The flow through the four cardiac valves will be investigated both in systole and diastole, resulting in the following eight measurements:

• Systolic and diastolic flow through the tricuspid valve
• Systolic and diastolic flow through the mitral valve
• Systolic and diastolic flow through the aortic valve
• Systolic and diastolic flow through the pulmonary valve

All the values are then integrated and a more reliable hemodynamic analysis is derived.

Broadly speaking, the blood flows could be analysed by considering the two following ratios: the atrio-ventricular ratio defined as

and the ventricular-arterial ratio:

Common clinical situations can therefore be described by comparing these two ratios, as presented in the following table:

Indeed, the comparison of these measurements with the unity may identify typical pathological features, helping the establishment of the diagnostic. Each one of the above-mentioned clinical situations can be interpreted as follows:

1. The amount of blood that is sent to the pulmonary circulation through the tricuspid and the pulmonary valve is similar to the amount of blood sent to the systemic circulation through the mitral and the aortic valve.


2. If R_TMF > 1 ±5% (and so probably R_PAF > 1 ± 5%), then the flow across the tricuspid valve is significantly higher than the flow across the mitral valve, which means that the right ventricle receives more blood than the left one. If the atria receives the same amount of blood we can deduce then that there is a left to right shunt. This parameter is crucial when evaluating atrial septum defects for example.


3. If R_TMF < 1 ±5% (and so probably R_PAF < 1 ± 5%), the phenomenon is inverted and in this case we can deduce that the shunt is right to left.


4. R_TMF is normal but R_PAF > 1 ± 5%. In this case there is more blood that passes through the pulmonary valve than through the aortic valve. One can deduce thus that there is a pulmonary regurgitation. This situation is found in post-operative patients who suffered from tetralogy of Fallot for instance.


5. R_TMF is normal but R_PAF < 1 ± 5%. This case is opposite to the previous situation, there is an aortic regurgitation. This phenomenon can be seen for instance in some hypertrophied cardiomyopathies where the pathological and thick septum shifts the aortic valve which, as a result, will not close accordingly.

Finally, R_PAF ratio may also help the expert to detect shunts ant to assess them as:

• trivial: if R_PAF is greater than 1.05 but lower than 1.5,
• mild: if R_PAF is greater than 1.5 but lower than 2,
• significant: if R_PAF is greater than 2 but lower than 3,
• severe: if R_PAF is greater than 3.

Second step: volume and contractility comparisons

Though comparing the flow across the various cardiac valves allows the cardiologists to detect some abnormal functions of the heart, other parameters may be useful to diagnose diseases that do not affect the blood circulation directly. In those cases, volumetric parameters may highlight atypical features of the myocardium and help the expert to identify cardiomyopathies or severe right-ventricle overload. For example,

• if the flows are normal and the volume is significantly large (+2 standard deviations), then the contractility of the muscle is abnormal. Actually, ventricles with poor contractility need to receive more blood in order to eject the same quantity of blood in the circulation system,
• if the flows are abnormal and the volume is significantly large, then a possible cause of the dilatation of the ventricle may be an abnormal blood circulation. This is typically the case in right-ventricular overload and in tetralogy of Fallot where the ventricle is dilated because of the increased pre-load. In atrial septal defect for example the abnormal flow is at the level of the atrio-ventricular valves whereas in post-operative tetralogy of Fallot the regurgitations of the pulmonary valve lead to a dilatation of the ventricle.

5. Conclusion

The main idea of knowledge-based models is to rely exclusively on standardised measurements that are acquired during clinical routine and to organise them into an in-silico algorithm developed from the methods that are used by the cardiologists.

A selection of relevant measurements that best describe the heart anatomy in children is presented in this document. The aim is to select and define wise clinical parameters that allow objective discrimination between the various pathologies. So far, cavity volumes, contractilities and flow measurements have been selected but other measurements might be added.

We have then proposed a preliminary knowledge-based model that allows the identification of typical clinical situations by analysing the patient's parameters. Moreover, such a method may also be applied to assess the severity of a disease by comparing the parameters to a standard scale.

However, the model might be improved by adding for instance more discriminant measurements to identify subtle pathological behaviours. This drives us to consider interesting challenges, such as the adjustment of the available imaging techniques to the specific case of paediatric diseases and to the specificities of the evaluation of the right-ventricle. Similarly, the decision tests may be refined in order to take into account all the possible clinical cases, from the most evident to the most subtle.

Disease-Based Models of Cardiac Electromechanics

1. Electromechanical model of the heart

1.1 Introduction

Several computational models of the heart have been proposed in the literature to reproduce the cardiac activity in a realistic way. Unfortunately, they tend to have an important number of parameters that are difficult to set. The purpose of the disease-based models is to simulate paediatric diseases by using a realistic computational model but whose parameters can be set from clinical and biological information. We decided therefore to rely on the approach proposed by () and developed at INRIA. Stated simply, it has a limited number of parameters, all related to clinically accessible measurements.

1.2 Anatomical model of the normal heart

The first step towards an in-silico heart model is to generate a synthetic anatomical representation of the normal organ. In this exploratory experiment we made use of the approach developed by (Sermesant et al., 2006), approach that we slightly adjusted to get a synthetic heart whose clinical parameters are close to the real clinical values.

The heart is modelled as a volumetric 3-D mesh based on tetrahedra. Contrary to the surface representations, the volumetric modelling allows the simulation of structural biomechanical phenomena, such as the muscle fibre contraction for example. Moreover, anatomical information like the variation of the cardiac muscle fibres across the myocardium for instance are naturally implemented by volumetric models. At last, clinical measurements are often preformed by considering the thickness of the cardiac muscle, measurements that can be easily simulated in this framework.

In this way, the creation of the anatomical model consists in four main steps:

• the design of the heart geometry,
• the synthesis of the cardiac muscle fibres related to the heart geometry,
• the segmentation of the resulting anatomy into regions of interest (according to the American Heart Association recommendations, (Cerqueira et al., 2002)),
• the creation of the 3-D mesh.

Several clinical studies (Mercier et al., 1982) have shown that the left ventricle (LV) can be well modelled by an ellipsoid truncated at the basal plane. On the other hand, the right ventricle presents with less regular shape and therefore there is no common agreement about a simplified geometrical representation (see next figure).

Short axis view of a normal heart (image acquired at Gaslini).
The right ventricle is at the left side of the images, while the left ventricle is at the right side.

In (Sermesant et al., 2006), the right ventricle (RV) is also represented by an ellipsoid. However, although this section of the heart has an ellipsoidal shape along the long axis, it looks like more "squared" than the left ventricle when seen in short axis. We therefore model the right ventricle by using a truncated super-ellipsoid adapted to our purposes. The super-ellipsoids extend the ellipsoid function and allow the creation of more rectangular geometries. In our application, they are defined by the following equations:

Formula of the ellipsoid (left) and of the super-ellipsoid (right) that models the left and right ventricles.

In the equation of the super-ellipsoid (the one to the right), the parameter n defines the degree of the XY-plane elements of the 3-D surface. The greater n is, the more rectangular the simulated shape is. Observe that the degree of the super-ellipsoid along the z-axis has a fixed value of 2 to constrain the simulated anatomical model to have an ellipsoidal shape along the base-to-apex direction.

2-D top view of the ellipsoid and super-ellipsoid models.
In red, an ellipsoid (d=2), in blue, super-ellipsoids whose degree is lower than 2 (d=0.5, 1 and 1.5),
in black, super-ellipsoids whose degree is higher than 2 (from d=2 to d=4, step of 0.5).

Two simulated heart anatomies.
Left panel: right ventricle is modelled by using a super-ellipsoid of degree d=2.5.
Right panel: right ventricle is modelled by using an ellipsoid.

By using this strategy, clinical parameters computed from the simulated heart are more realistic. For example, the ratio between the right-ventricular volume and the left-ventricular volume is equal to 0.67 when the bi-ellipsoidal model is used, whereas it is equal to 0.90 when the right-ventricle is represented by a super-ellipsoid with a degree of 2.5. Observed volume ratios are bounded by 0.9 and 1.1 in normal subjects.

It is worth mentioning that, as in (Sermesant et al., 2006), this anatomical model can be easily adapted to pathological cases by modifying global parameters such as the radii of the cavities, the thickness of the myocardium, the length of the heart, etc.

The second step consists in generating the muscle fibres tailored to the synthetic myocardium created during the previous stage. As stated in (Ubbink et al. 2006), the organisation of the myocytes (the muscle fibres of the heart) and their orientation play a crucial role in the function of the normal heart. Contractility anomalies for example often raise from pathologies that affect the myocyte cells.

It is thus essential to correctly model the fibres and their variation across the cardiac wall to get reliable simulations. Fortunately, many dissection analyses and, more recently, diffusion tensor imaging have shown that the fibre orientation is quite constant among the individuals of a same species. In humans, the myocytes are organised according to a common pattern:

• the angle with the short axis varies across the cardiac wall from around +70° on the endocardium to around -70° on the epicardium (Hsu and Henriquez 2001, Guccione and McCulloch, 1991),
• at the mid-wall they are almost horizontal in the short axis plane.

Simulated orientations of the cardiac muscle fibres.

The third step of the algorithm aims at assigning at each tetrahedron of the volumetric anatomical model a label corresponding to a standard anatomical region. The main rationale here is that accurate calibration, estimation and analysis of the model is made easier by subdividing it into different anatomical regions. Moreover, clinical parameters are often measured in selected region of interests since the pathologies may be focal and affect a localised part of the heart only.

In this way, the synthetic anatomical model is divided into the segments proposed by the American Heart Association (Cerqueira et al., 2002). This subdivision will allow the simulation of local abnormalities in pathological cases, such as the focal thickening of the septum and local fibre disarray in hypertrophied cardiomyopathy.

AHA segmental division of the simulated anatomy.

Finally, the volumetric mesh is created. The marching cube algorithm (Lorensen and Cline, 1987) is applied to extract a triangulated surface from the anatomical image created during the first step. The result is then used to create the tetrahedral volumetric mesh (INRIA GHS3D software). Finally, fibre orientations and anatomical subdivisions are assigned to each tetrahedron of the mesh through rasterisation (Sermesant et al., 2003).

Observe that the simulated geometry corresponds to the passive state of the heart, where the fibres are neither stretched by any pre-load nor contracted.

INTERACTIVE 3D MESH

Synthetic anatomical model.
Left panel: the modelled anatomy with simulated AHA segments and muscle fibres.
Right panel: interactive 3-D mesh (left mouse button: rotation, right mouse button: zoom).

1.3 Transmembrane potential simulation

Many different models for the simulation of the electrophysiology of the heart have been proposed. Two main categories can be identified: the biophysical or ionic models, where the heart is described at the cellular level, and the phenomenological models, where the cardiac system is expressed at a more macroscopic level. Contrary to the biophysical models, the phenomenological methodology does not compute the ion concentrations that govern the electrophysiology of the heart but aims at solving equations that describe global phenomena. The resulting systems are thus simpler and, above all, depend on parameters that can be measured during clinical evaluation.

In our approach, the method proposed in (Sermesant et al., 2006) is applied. It consists of a simplified phenomenological model based on the version proposed by (Aliev and Panfilov, 1996) of the FitzHugh-Nagumo equations. This model is well adapted to the dynamics of the cardiac electrical potential and captures the main biological features of the heart function, which are:

• the cell is activated only for a stimulus larger than a certain threshold,
• the shape of the transmembrane potential does not depend on the stimulus,
• the cell cannot be excited during a refractory period,
• any cell can be stimulated.

(Sermesant et al., 2006) modified the Aliev and Panfilov model by ignoring the pacing frequency of the transmembrane potential duration, which gave the following equations:

The orientation of the fibres is modelled as an anisotropic 3x3 conductivity tensor D in an orthonormal basis. The first eigenvector is along the local fibre orientation f, the conductivity in the fibre direction is expressed by the constant d₀ and r is the conductivity anisotropy ratio in the transverse plane, typically four times as low as d₀:

Top: simulated transmembrane potential. Time is normalised to obtain a duration equal to 1 time unit.
Bottom: measured transmembrane potential on a frog cardiac cell.

Finally, these equations are integrated into the 3D volumetric mesh. A normalised transmembrane potential of 1.0 is imposed at the nodes of the endocardium surfaces as an initial condition. The resulting propagation can then be visualised with the isochrones and the different polarisation times are represented by colours.

1.4 Biomechanics simulation

In a general point of view, the myocardium is made up of active nonlinear anisotropic visco-elastic materials. Its constitutive law is complex and must include an active element, which simulates the contraction of the muscle and is controlled by the transmembrane potential, and a passive element representing the mechanical elasticity.

Several constitutive laws have been proposed in the literature (Hunter and Smaill, 1988; Humphrey et al., 1990; Hunter et al., 1997; Nash, 1998; Häfner et al., 2002) but as the Health-e-Child project aims at developing advanced model integrating vertical data (from genetics to physiological), the hierarchical approach introduced in (Bestel et al., 2001) and adapted by (Sermesant et al., 2006) will be used.

In (Bestel et al. 2001), the authors proposed to model the contraction from the nano-motors scale and to build up the macroscopic constitutive law by representing the phenomena encountered at different scales. The model is based on the Hill-Maxwell rheological scheme: the muscle is represented by the serial combination of a contractile element with an elastic element for the isovolumetric contraction and by a parallel element mainly representing the elasticity of the material.

Scheme of the rheological model. E_c: active contractile element, E_p: passive element.

• an active contractile element, E_c, which creates the stress tensor σ_c and is controlled by the normalised transmembrane potential u.
• a passive parallel element, E_p, which is anisotropic linear visco-elastic and creates a stress tensor σ_p

The electromechanical coupling is performed by using a simple ordinary differential equation. This is motivated by the fact that few clinical parameters are available and the model must thus be as simple as possible. In this way, we choose to compute the contraction intensity from the normalised transmembrane potential u by solving:

Moreover, because the normalised transmembrane potential is bounded by 0 and 1 and the changes on depolarisation and repolarisation are abrupt, u can be approximated by the values 0 and 1, using the current computed value thresholded at 0.5. The differential equation is then directly solvable and the contractility of the heart is controlled by the following coupling equation system:

T_d is the depolarisation time, T_r the repolarisation time, HP the heart period, α_c the contraction rate, α_r the relaxation rate, σ₀ is the maximum contraction and σ_r=σ_c(T_r). α_c and α_r control the contraction stress increase and decrease while T_d and T_r allow the simulation of the delay between the electrical signal and the biomechanical answer.

The 3-D anisotropic contraction stress tensor that takes into account the fibre orientation is obtained from the coupling equations. The force vector F_c is then estimated and integrated in the dynamics equation:

In this equation U is the displacement vector, M the diagonal mass matrix (mass lumping), C a diagonal damping matrix, F_b the different external loads from the boundary conditions and K the anisotropic linear elastic stiffness matrix. Since the material under consideration is linear elastic but anisotropic, K is constant in small displacement formulation. Its construction is based on the finite element method with linear tetrahedral elements, with the derivation of displacements U into the linearised strain tensor and with the Hookean constitutive law between Cauchy stress tensor and the strain tensor.

1.5 Boundary conditions

We aim at simulating the cardiac function during one or several cardiac cycles. The interaction of the myocardium with the blood is therefore crucial and the strategy described in (Sermesant et al., 2006) is implemented and slightly modified for our purposes. The four different phases of the cardiac cycle are modelled and different boundary conditions are applied at each phase:

• filling: during this phase the ventricles are filled by the blood that comes from the atria. We apply thus to the vertices of the endocardia a pressure whose intensity is equal to the mean pressure of the atria. In order to simulate the atrial contraction, this pressure is augmented during the P wave.
• isovolumetric contraction: here the myocardium begins to contract but the pressure inside the ventricles is not sufficiently high to open the arterial valves: the blood volume is thus constant. To simulate this phase the approach presented in (Sermesant et al., 2006) is used: a penalty constraint is applied to the vertices of the endocardia to keep the ventricle volume constant in spite of the contraction. This constraint is estimated at each time step and its value corresponds to the ventricular pressure during this cardiac phase.
• ejection: once the pressure inside the ventricles is high enough to open the arterial valves, the blood flows to the arteries. The pressure applied to the vertices of the endocardia is therefore equal to the mean pressure of the aorta (for the left ventricle) and the pulmonary artery (for the right ventricle).
• isovolumetric relaxation: the myocardium begins to relax but the pressure is too high and the atrio-ventricular valves are still closed: the blood volume is constant. The same strategy as the one used during the isovolumetric contraction is applied.

In the current implementation, atrial pressures have two values (baseline and atrial contraction) and the arterial pressures have a constant value. The transition algorithm presented in (Sermesant et al., 2006) has been adapted to support pathological behaviours and is performed as follows:

• from filling to isovolumetric contraction: when the contraction starts, the blood flow changes sign. The atrio-ventricular valves close and the isovolumetric contraction starts. However, this transition is performed only if the cardiac myocytes are being activated by the transmembrane potential, that is when the heart is contracting itself.
• from isovolumetric contraction to ejection: the arterial valves opens as soon as the ventricular pressure is higher than the arterial one
• from ejection to isovolumetric relaxation: when the relaxation starts, the blood flow changes sign. The valves close and the next isovolumetric relaxation starts. However, this transition is performed only if the cardiac myocytes are coming back to their original state, that is when the heart is relaxing itself.
• from isovolumetric relaxation to filling: when the ventricular pressure is lower than the atrial one, the atrio-ventricular valve opens and the filling starts.

Next sections present in more details the results of first simulations based on this model of the normal heart and paediatric diseases.

1.6 References

• Aliev, R., Panfilov, A., 1996. A simple two-variable model of cardiac excitation. Chaos, Solitons, Fractals, vol. 7, no. 3, pp. 293-301.
• Bestel, J., Clement, F., Sorine, M., 2001. A biomechanical model of muscle contraction. Lecture Notes in Computer Science. Berlin, Germany: Springer-Verlag, vol. 2208, Medical Image Computing and Computer-Assisted intervention (MICCAI 01), pp. 1159-1161.
• Cerqueira, M. D., Weissman, N. J., Dilsizian, V., Jacobs, A. K., Kaul, S., Laskey, W. K., Pennell, D. J., Rumberger, J. A., Ryan, T., Verani, M. S., 2002. Standardized myocardial segmentation and nomenclature for tomographic imaging of the heart: a statement for healthcare professionals from the Cardiac Imaging Committee of the Council on Clinical Cardiology of the American Heart Association. Circulation, Vol. 105, No. 4., pp. 539-542.
• Guccione, J., McCulloch, A., 1991. Theory of Heart: Biomechanics, Biophysics, and Nonlinear Dynamics of Cardiac Function. Springer-Verlag, pp. 121-144 (Ch. Finite element modeling of ventricular mechanics)
• Häfner, J., Sachse, F., Sansour, C., Seemann, G., Dössel, O., 2002. Hyperelastic description of elastomechanic properties of the heart: a new material law and its application. Biomedizinische Technik, vol. 47-1/2, pp. 770-773.
• Hsu, E., Henriquez, C., 2001. Myocardial fiber orientation mapping using reduced encoding diffusion tensor imaging. Journal of Cardiovascular Magnetic Resonance 3, 325-333.
• Humphrey, J., Strumpf, R., Yin, F., 1990. Determination of a constitutive relation for passive myocardium: I. A new functional form. ASME J. Biomechan. Eng., vol. 112, pp. 333-339.
• Hunter, P., Smaill, B., 1988. The analysis of cardiac function: a continuum approach. Progr. Biophys. Mol. Biol., vol. 52, pp. 101-164
• Hunter, P., Nash, M., Sands, G., 1997. Computational electromechanics of the heart. Computational Biology of the Heart. New York: Wiley, ch. 12, pp. 345-407
• Lorensen, W., Cline, H., 1987. Marching cubes: a high resolution 3D surface reconstruction algorithm. Computer Graphics (Proc. SIGGRAPH), vol. 21, no. 4, pp. 163-169.
• Mercier J., Disessa T., Jarmakani J., Nakanishi T., Hiraishi S., Isabel-Jones J., Friedman W., 1982. Two-dimensional echocardiographic assessment of left ventricular volumes and ejection fraction in children. Circulation, Vol. 65, No. 5., pp. 962-969.
• Nash, N., 1998. Mechanics and material properties of the heart using an anatomically accurate mathematical model. Ph.D. dissertation, Univ. Auckland, Auckland, New Zeland.
• Sermesant, M., Forest, C., Pennec, X., Delingette, H., Ayache, N., 2003. Deformable biomechanical models: Application to 4D cardiac image analysis. Med. Image Anal. 7, 475-488.
• Sermesant, M., Delingette, H., Ayache, N., 2006. An Electromechanical Model of the Heart for Image Analysis and Simulation. IEEE Trans. Med. Imag., vol. 25, no. 5, pp. 612-625.
• Ubbink, D., Bovendeerd, P., Delhaas, T., Arts, T., van de Vosse, F., 2006. Towards model-based analysis of cardiac MR tagging data: relationn between left ventricular shear strain and myofiber orientation. Medical Image Analysis. Special Issue on FIMH.

Disease-Based Models of Cardiac Electromechanics

2. Simulation of the normal heart

The normal heart is simulated first in order to validate the model and to calibrate the parameters to the real values. The resulting simulation protocol will then be considered as reference when tailoring the algorithm to the pathologies under study.

2.1 Parameters of the normal heart

To simulate the normal heart, the geometrical parameters are defined as follows:

• The degree of the super-ellipsoid that models the right-ventricle is chosen equal to 2.5 in order to have a simulated anatomy consistent with the clinical observations (see next table). The right ventricle is therefore defined by the equation:

• The orientation of the cardiac muscle fibres varies from +90° on the endocardium to -90° on the epicardium (Sermesant et al., 2006)

INTERACTIVE 3D MESH

Geometrical model of the normal heart.
Left panel: simulated geometry, the black lines represent the cardiac muscle fibres.
Right panel: interactive 3-D mesh, left click: rotation, right click: zoom.

The electrical parameters are defined according to (Sermesant et al., 2006) and to (Aliev and Panfilov, 1996): e = 0.01, k = 8, a = 0.15 (more information). The actual transmembrane potential E in mV is obtained by using the conversion formula E = 100 * u - 80 and the actual time is computed by scaling the normalised time used in the electro-physiological equations (where the duration of the simulated action potential is equal to 1.0) to get a more realistic value of 300 ms.

Finally, the conductivity d₀ in the fibre direction is equal to 1.0 and the conductivity in the transverse plane r is equal to 0.25.

The biomechanical parameters are tailored to the new geometry and are chosen in order to simulate a normal cardiac function in children. In this way, the contraction rate α_c is now equal to 30, the relaxation rate α_r to 20 and the maximum contraction σ₀ to 0.025. A delay between the electrical signal and the biomechanical answer is also considered: T_d = 0.01 and T_r = 0.3.

In this model the pressure in the atria and in the arteries are constant all along the cardiac cycle. We therefore set these parameters as close as possible to their mean values during the cardiac cycle. Nevertheless, a compromise between numerical stability and respect of the actual values had to be found. Thus, pressure in the left atrium is set equal to 15 mmHg, in the aorta to 60 mmHg, in the right atrium to 7.5 mmHg and in the pulmonary artery to 30 mmHg.

The following table provides a unified view of the parameters of the simulation.

2.2 Results and discussion

Four cardiac cycles are simulated but only the last three are considered for the evaluation. The first cycle is a transitional state of the system and may thus bias the results.

At each cardiac cycle and for each ventricle we calculate four clinical measurements:

• end-diastolic volume (EDV): volume of the cavity when the heart quits the filling phase and starts the isovolumetric contraction,
• end-systolic volume (ESV): volume of the cavity when the heart quits the ejection phase and starts the isovolumetric relaxation,
• stroke volume (SV): amount of blood ejected during a cardiac cycle, defined by SV = EDV - ESV,
• ejection fraction (EF): indicator of global contraction, defined by EF = SV / EDV.

Since the system converges after two cardiac cycles, only the fourth and last simulated cycle is reported in next table. The diagram of the volume variations over the entire simulation is shown in the figure below.

Simulated volume variations of the normal heart. In red, left ventricle. In blue, right ventricle.

Left ventricle. According to (Helbing et al., 1995), the simulated end-diastolic and end-systolic volumes of the left ventricle are consistent with clinical observations in healthy children (LVEDV = 89 ± 26 mL, RVEDV = 26 ± 9 mL). Left ejection fraction can be seen as normal or slightly subnormal. In (Helbing et al., 1995), LVEF = 70% ± 6.

Right ventricle. If the simulated end-diastolic volume is consistent with clinical observations (RVEDV = 92 ± 25 mL according to (Helbing et al., 1995)), the simulated end-systolic volume is moderately high and cannot be considered as a normal value in children (39.77 mL versus 27 ± 9 mL (Helbing et al., 1995)). Consequently, the simulated ejection fraction is lower than the observed ejection fraction in healthy children (47.26% versus 70% ± 4, (Helbing et al., 1995)).

The ratio between the end-diastolic volumes of the right and left ventricles (equal to 0.9) is consistent with the observed values (between 0.9 and 1.1 depending on the studies (Graham et al., 1973; Helbing et LL., 1995; Lorenz 2000)).

Finally, the simulation has been assessed visually. Both ventricles stay synchronised during all the cardiac phases, and no abnormal septal movement is apparent. The simulated normal heart is shown in the following video.

2.3 Conclusion and future works

The disease-based model previously presented has been calibrated by simulating the normal heart. The parameters of the simulation have been set according to clinical information and, for those where no clinical data was available, to get simulated clinical parameters consistent with actual measurements.

The resulting simulations are promising. The values of the simulated cardiac parameters related to the left ventricle fall in the range of the observed values in healthy children, whereas those related to the right ventricle are slightly different but close to the clinical references.

The inconsistencies in the right-ventricle parameters are mainly due to a limitation of the anatomical model. Indeed, to build the synthetic anatomy, the geometrical primitives are cut at the basal plane. Part of the right ventricle, in particular the outflow tract, is therefore ignored by our 3-D model. The volume of the right cavity is thus under-estimated and the right-ventricular wall truncated.

A way to solve this issue would consist in modelling the heart anatomy by cutting the primitives by the valvar plane instead of the basal plane. The use of advanced techniques such as those presented in the last section of this document constitutes another solution currently under study in Health-e-Child.

Another limitation raised by this first simulation is the use of constant pressures as boundary conditions. Indeed, in the real heart these pressures vary significantly during the cardiac cycle. The reliability of the simulations would therefore be improved if such phenomena were implemented. Ejection fraction for example may be more accurate. To this aim, we plan to implement advanced techniques that simulate the pressure variations along the cardiac cycle, such as the Windkessel model (Sermesant et al., 2006b).

Yet, promising results have been obtained despite these limitations. They show that it is possible to simulate paediatric cardiac diseases by using 3-D electromechanical models of the heart. This drove us to perform preliminary simulations of right-ventricle overload and dilated and hypertrophic cardiomyopathy. The first results are presented in the following sections.

2.4 References

• Aliev, R., Panfilov, A., 1996. A simple two-variable model of cardiac excitation. Chaos, Solitons, Fractals, vol. 7, no. 3, pp. 293-301.
• Graham, T. P., Jarmakani, J. M., Atwood, G. F., Canent R. V., 1973. Right Ventricular Volume Determinations in Children: Normal Values and Observations with Volume or Pressure Overload. Circulation, Vol. 47, No. 1., pp. 144-153.
• Helbing, W. A., Rebergen, S. A., Maliepaard, C., Hansen, B., Ottenkamp, J., Reiber, J. H., de Roos A., 1995. Quantification of right ventricular function with magnetic resonance imaging in children with normal hearts and with congenital heart disease. Am Heart J, Vol. 130, No. 4., pp. 828-837.
• Lang, R. M., Bierig, M., Devereux, R. B., Flachskampf, F., A., Foster, E., Pellikka, P., A., Picard, M., H., Roman, M., J., Seward, J., Shanewise, J., S., Solomon, S., D., Spencer, K., T., Sutton, M., St J., Stewart, W., J., 2005, Recommendations for Chamber Quantification: A Report fromt the American Society of Echovardiography's Guidelings end Standards Committee and the Chamber Quantification Writing Group, Developed in Conjunction with the European Association of Echocardiography, a Brandch of the European Society of Cardiology. J. Am. Soc. Echocardiogr, no 18, pp. 1440-1463
• Lorenz, C. H., 2000. The Range of Normal Values of Cardiovascular Structures in Infants, Children, and Adolescents Measured by Magnetic Resonance Imaging. Pediatric Cardiology, Vol. 21, No. 1., pp. 37-46.
• Pfisterer, M. E., Battler, A., Zaret, B. L., 1985. ange of normal values for left and right ventricular ejection fraction at rest and during exercise assessed by radionuclide angiocardiography. Eur Heart J, vol. 6, no. 8, pp. 647-655
• Sermesant, M., Delingette, H., Ayache, N., 2006. An Electromechanical Model of the Heart for Image Analysis and Simulation. IEEE Trans. Med. Imag., vol. 25, no. 5, pp. 612-625.
• Sermesant, M., Moireau, P., Camara, O., Sainte-Marie, J., Andriantsimiavona, R., Cimrman, R., Hill, D. L. G., Chapelle, D., Razavi, R., 2006b. Cardiac function estimation from MRI using a heart model and data assimilation: Advances and difficulties. Medical Image Analysis, vol. 10, pp. 642-656.

Disease-Based Models of Cardiac Electromechanics

3. Simulation of right-ventricle overload

3.1 Cardiac parameters for right-ventricle overload

Diagram of atrial septal defect

Right-ventricle overload (RVO) is a feature of some cardiac diseases. It is characterised by an increased pre-load of the right ventricle due to either an abnormally high income of blood from the right atrium (in atrial septal defects for instance) or to pulmonary valve regurgitation (in post-operative patients that suffered from the tetralogy of Fallot for example). The increased pre-load leads to a dilatation of the right ventricle heart and may locally alter the muscle integrity. Also, boundary conditions at the valves are modified because of the abnormal blood flow.

In this chapter a simulation of RVO is presented, where only geometrical and boundary condition changes are considered. In particular, we aim at simulating RVO diseases that are characterised by increased atrio-ventricular pressures, as in atrial septal defects. Clinical observations are used to set the parameters of the 3-D model.

Clinically, one of the most relevant features of RVO is the deformation of the right ventricle. Due to the increased pre-load, the ventricle is dilated but without relevant increase of myocardium mass. The end-diastolic volume is therefore higher than in normal subjects (128 mL/m² with respect to 70 mL/m² in normal children (Graham et al., 1973)), just like the volume ratio between right and left ventricle (2.3 versus 1 (Graham et al., 1973)).

The dilatation of the right ventricle is simulated by a higher major axis of the super-ellipsoid that models this cavity (+29%). Moreover, in order to simulate the stretching of the myocardium, the synthetic fibre directions are made more horizontal, varying from +60° on the endocardium to -60° on the epicardium. Finally, the conservation of the mass is simulated by thinning the right-ventricular myocardium (-20%). The mass of the simulated pathological geometry is then almost equal to the one of the synthetic normal heart (98%). The figures below show the resulting anatomy of RVO and a comparison with the geometrical model of the normal heart.

RVO geometrical model.
Left panel: simulated RVO anatomy. Black lines represent the muscle fibres.
Right panel: comparison between RVO (in red) and normal (in blue) heart geometry. The RVO right cavity is larger than the normal cavity, with thinner myocardium.

Boundary conditions of the right ventricle are also altered. In particular, the filling pressure during the filling phase is higher than in normal subjects. In this way, the pressure of the right atrium is set equal to 12.8 mmHg, which is greater than the normal pressure (7.5 mmHg) but slightly lower than the pressure of the left atrium (15 mmHg).

Next table summarises the parameters used to simulate RVO. Normal parameters are reported for comparisons, the differences are highlighted in blue.

The other parameters, in particular the electromechanical parameters, are kept unchanged. (see normal values).

3.2 Results and discussion

Four cardiac cycles have been simulated. As for the simulation of the normal heart, the steady state of the system is reached at the second cycle. We therefore report the quantitative results computed from the latest cycle only. The following figure illustrates the volume variations of the simulated RVO heart.

Simulated volume variations. Left panel: simulation of right-ventricle overload.
Right panel: simulation of a normal heart. The scales are aligned for comparison.

These results are compared with normal values and with clinical observations in children.

In addition to the dilatation of the right ventricle, RVO is characterised by an increased RVEF/LVEF ratio, but with normal right-ventricle ejection fraction (RVEF) (Graham et al., 1973), and increased end-systolic and end-diastolic volumes (Helbing et al., 1995).

From the simulations RVEF/LVEF ratio amounts to 0.88 for RVO and to 0.77 for the normal heart, that is an increase of 14%. This increment is in agreement with clinical observations: RVEF/LVEF is equal to 0.99 in healthy children and to 1.17 in RVO subjects, that is an increment of 17%.

On the other hand, the right ejection fraction has not significantly increased with respect to the normal heart (53.79% in contrast with 48.76%). This progression is consistent with the clinical surveys. However, the simulated ejection fraction is still low in comparison with ground-truth measurements (54% versus 70%), as in the simulation of the normal heart.

Finally, the simulated end-diastolic and end-systolic volumes of the right ventricle are significantly higher than those of the simulated normal heart, with a volume ratio equal to 1.9 at end diastole and to 2.3 at end systole (0.9 and 1.3 in the simulated normal heart). These results are in line with clinical measurements.

3.3 Conclusion

This experiment constitutes an attempt to simulate RVO. By slightly changing the heart geometry and by modifying the boundary conditions at the level of the tricuspid valve we managed to retrieve clinical features of this disease.

Clinical cardiac parameters have changed accordingly with respect to real observations: the ejection fraction of the right ventricle was constant whereas the ratio RVEF/LVEF and the right-ventricle volumes at end diastole and end systole increased. Though the dilatation of the right ventricle has been simulated a priori during the modelling of the pathological heart anatomy, the increased pre-load amplified the phenomenon (RV volume / LV volume ratio increases from 1.58 to 1.92). This result is promising since it demonstrates the capabilities of the 3-D model to adapt himself to abnormal conditions.

Furthermore, clinical studies have shown that the contractility properties of the myocardium vary according to the filling pressure. Owing to Laplace's law, the more the muscle fibres are stretched, the higher the contractility is. The ejection fraction increases, till a certain limit, with the filling pressure. In our simulation, this property has not been implemented explicitly but, since the synthetic cardiac muscle are modelled by using visco-elastic methods, the increased stretching of the finite elements have lead to a major contractility which results in an augmentation of the ejection fraction (from 48.8% to 53.8%). However, this phenomenon is completely passive. No active biological phenomenon has been considered in this experiment.

However, there is still some limitations. As in the simulation of the normal heart, the ejection fraction of the right ventricle is too weak to be considered as normal in children. Further analyses of the right-ventricular geometry and data assimilation procedures are required to improve the reliability of the simulation.

Finally, other right-ventricle overload diseases are planed to be simulated, such as the post-operative tetralogy of Fallot or partial anomalous pulmonary venous return.

3.4 References

• Graham, T. P., Jarmakani, J. M., Atwood, G. F., Canent R. V., 1973. Right Ventricular Volume Determinations in Children: Normal Values and Observations with Volume or Pressure Overload. Circulation, Vol. 47, No. 1., pp. 144-153.
• Helbing, W. A., Rebergen, S. A., Maliepaard, C., Hansen, B., Ottenkamp, J., Reiber, J. H., de Roos A., 1995. Quantification of right ventricular function with magnetic resonance imaging in children with normal hearts and with congenital heart disease. Am Heart J, Vol. 130, No. 4., pp. 828-837.
• Lorenz, C. H., 2000. The Range of Normal Values of Cardiovascular Structures in Infants, Children, and Adolescents Measured by Magnetic Resonance Imaging. Pediatric Cardiology, Vol. 21, No. 1., pp. 37-46.
• Pfisterer, M. E., Battler, A., Zaret, B. L., 1985. Range of normal values for left and right ventricular ejection fraction at rest and during exercise assessed by radionuclide angiocardiography. Eur Heart J, vol. 6, no. 8, pp. 647-655

Disease-Based Models of Cardiac Electromechanics

4. Simulation of dilated cardiomyopathy

4.1 Cardiac parameters for dilated cardiomyopathy

Dilated cardiomyopathy (DCM) is characterised by left-ventricular or, to a lesser extent, bi-ventricular dilatation with severely depressed systolic function due to impaired myocardial contractility (Dec and Fuster, 1994). At a microscopic level, hypertrophy and degeneration of myocytes with varying degree of interstitial fibrosis are common features. In some cases, small clusters of lymphocytes are also present (Dec and Fuster, 1994).

The pathogenesis is unclear and DCM can be the result of any variety of cytotoxic, metabolic, immunological, genetic or infectious mechanisms (Dec and Fuster, 1994) that damage the myocardium. The main clinical symptom of DCM is left-sided heart anomaly or indeed failure. Premature ventricular beats, arrhythmias and decompensated heart failure are the most relevant physical findings.

Unfortunately, there are few well-defined features of DCM. The most relevant ones are the dilatation of the heart and an abnormal systolic function but they are very heterogeneous and their severity varies significantly among the patients. We will however rely on these attributes to perform a simulation of DCM.

To simulate DCM, the normal heart anatomy must be modified so as to enlarge the ventricles. In this experiment, both ventricles are dilated. Left ventricle is made more spherical (Dec and Fuster, 1994) by increasing the minor axis of the ellipsoid modelling this cavity (+45%). Similarly, the right ventricle is dilated by increasing its minor axis (+28%) and its depth: the right-ventricular apex is at the same depth as the left-ventricular apex.

The volumes of the dilated cavities are about twice as great as those of the normal heart (see next table).

Observe that in DCM, contrary to right-ventricle overload, the myocardium is dilated but its thickness is normal. The pathological heart is therefore heavier than the normal organ. In our simulation, the thickness of the cardiac wall is normal but because of the dilatation its weight amounts to 1.4 times the weight of the normal heart.

Finally, fibre orientations are assumed to be normal, varying from +90° on the endocardium to -90° on the epicardium.

Synthetic anatomy of DCM. Black lines correspond to the cardiac fibre directions.

Comparison between DCM (in red) and normal (in blue) anatomy. Both ventricles of the DCM heart are dilated, and right-ventricular and left-ventricular apices are at same depth.

From a functional point of view, DCM is characterised by an impaired contractility caused by the degeneration of the myocytes and other microscopic abnormalities. It is therefore crucial to integrate this feature into the 3-D model.

As the 3-D model of the heart is controlled by macroscopic parameters only, the microscopic abnormalities are modelled by reducing the maximum contraction parameter &sigma₀ (more information here), from &sigma₀ = 0.025 to &sigma₀ = 0.01.

No more parameters are modified for this simulation.

Arteries and atria are supposed to be healthy: boundary conditions stay unchanged.

The table below summarises the parameters of the simulation. Normal values are mentioned for comparisons. The parameters that have been modified are highlighted in blue.

The other biomechanical parameters are kept unchanged (see normal values).

4.2 Results and discussion

For the normal heart and for right-ventricle overload, four cardiac cycles are simulated. The quantitative measurements estimated from the last simulated cycle are reported in the next table. The following figure shows the volume variations along all the simulated cycles.

Simulated volume variations. Left panel: simulation of dilated cardiomyopathy.
Right panel: simulation of a normal heart. The scales are aligned for comparison.

Since the abnormal geometry has been explicitly modelled, the only well-known clinical DCM feature that can come up from the simulation is the low ejection fraction. The simulated DCM heart has a LV ejection fraction equal to 22.69 % and a RV ejection fraction equal to 15.90 %, which are values likely to be found in patients with severely impaired systolic function (Dec and Fuster, 1994).

4.3 Conclusion

A simulation of dilated cardiomyopathy has been presented. The heart geometry has been modified according to clinical observations and the impaired contractility has been modelled by reducing the maximum contraction parameter.

The results are promising since they correspond to clinical observations: DCM leads to a significant reduction of the ejection fraction of the affected ventricles (Dec and Fuster, 1994).

However, by analysing the simulated volume variations, an abnormal phenomenon can be identified at the beginning of the filling and the ejection phases of the right ventricle. In next figure we have magnified the diagram of the volume variations to see a single cycle only, and the transitions are shown through the vertical bars. We can see in that way that when the filling starts the volume of the right ventricle first decreases and then increases.

Volume variations over one cardiac cycle. The dashed lines represent the phase transitions of the right-ventricle. The diagram starts at the filling phase, then comes the iso-volumetric contraction, then the ejection, followed by the iso-volumetric relaxation, and so on.

The visual analysis of the simulation sequence provides a possible reason to this issue: the filling of the dilated left ventricle makes the intra-ventricular septum move towards the right ventricle. However, the relaxation rate of the right ventricle is not sufficient enough to counter-balance the movement of the septum and the volume of the cavity decreases briefly as if there is a regurgitation (see video below).

These anomalous volume variations at the beginning of the ejection and filling phases are due to limitations of our heart model. In particular, regurgitations are not implemented and the phase transitions are performed in a sequential way. The abnormal closure of the valves during relaxation for example, which corresponds to a brief iso-volumetric relaxation, cannot be simulated with our current implementation where iso-volumetric relaxation can take place only after ejection and before filling. Impoving the DCM disease-based models thus requires a better modelling of the valves and the regurgitations.

4.4 References

• Dec, G. W. and Fuster, V., 1994. Idiopathic Dilated Cardiomyopathy. N Engl J Med, Vol. 331, No. 23., pp. 1564-1575.

Disease-Based Models of Cardiac Electromechanics

5. Simulation of hypertrophic cardiomyopathy

5.1 Cardiac parameters for hypertrophic cardiomyopathy

Like dilated cardiomyopathy (DCM), idiopathic hypertrophic cardiomyopathy (HCM) is a disease that affects the correct function of the myocardium. However, contrary to DCM, HCM is characterised by a focal or global thickening of the myocardium without enlargement of the cavities and, at a microscopic level, by myocardial fibre disarray and fibrosis.

HCM may be asymptomatic and diagnosed only after a genetic survey performed when a family member of the patient suffers from this disease. When HCM is symptomatic, the main clinical features that characterise it are congestive heart failures and heart rhythm disturbances which may lead to sudden death.

The pathogenesis is unknown but genetic mutations are more and more identified as possible reasons of HCM (Wigle et al., 1995, Maron, 2002). That is why genetic surveys are often realised when a patient is found to suffer from HCM. However, other factors may cause HCM. Indeed, they may be the result of congenital metabolic diseases (mitochondrial disorders) or of overload pathologies like for instance the Pompe disease (a glycogen storage disease) or Friedreich's ataxia (an inherited neuro-muscular disease), both affecting the myocardium. HCM may also be found in syndromes such as the Noonan syndrome, a congenital genetic condition.

There exist various types of HCM, with or without outflow obstruction. Since about 75% of the patients do not have any valvar obstruction (Maron, 2002), we chose to simulate this category of HCM first. Moreover, global HCM without ventricular dilatation is considered here, where the left ventricle is entirely affected by the disease but the right ventricle is healthy. In the following the parameters of the simulation are briefly presented.

Contrary to right-ventricle overload and dilated cardiomyopathy, there is no typical geometry of hypertrophic heart. The only common feature is a focal or global thickening of the myocardium. In this experiment, the normal geometry is modified by thickening the left ventricle myocardium (+28%), without cavity enlargement. The fibre orientations stay normal, i.e. from +90° to -90°, as well as the geometry of the right ventricle.

Geometrical model of HCM.
Left panel: simulated HCM anatomy. Black lines correspond to the cardiac fibre directions.
Right panel: Comparison between HCM (in transparent red) and normal (in blue) geometry.
The HCM left-ventricular myocardium is thicker than in the normal heart.

As stated in (Wigle et al., 1995), HCM is characterised at a microscopic level by a myocardial fibre disarray and fibrosis. This, along with the increased thickness of the cardiac muscle, leads to a stiffer muscle and an impaired cardiac relaxation affecting the diastolic function. In particular, relaxation rate and volume filling is reduced, which results in a slower myocardium relaxation.

To simulate these observations, the relaxation rate of the 3-D model is decreased (from α_r=20 in the normal heart to α_r=10) (more information here) and the maximum contraction is reduced from &sigma₀ = 0.025 to &sigma₀ = 0.022

Only these biomechanical parameters are explicitly modified since:

• stiffness will be implicitly increased owing to the thickening of the myocardium,
• electrical propagation speed will be implicitly reduced because of the enlargement of the myocardium of the left ventricle (the electrical waves must cross more tetrahedra due to the thickening of the cardiac muscle).

Since we aim at simulating hypertrophic cardiomyopathy without valve obstruction and with normal atria, the boundary conditions are kept at their normal values.

Next table summarises the parameters of the HCM heart. Normal values are mentioned for comparisons. The differences are highlighted in blue.

The other biomechanical parameters stay unchanged (see normal values).

5.2 Results and discussion

Four cardiac cycles have been simulated. The steady state of the system has been reached after the second cycle. We therefore report in the following table only the results of the last simulated cycle.

Simulated volume variations. Left panel: simulation of hypertrophic cardiomyopathy. Right panel: simulation of a normal heart.
The scales are aligned for comparison.

The main feature of the simulated volume variations is the slight increase of left-ventricular ejection fraction (from 63.2 % to 65.4 %, that is an increase of 3%), whereas the right-ventricular ejection fraction stays constant. This result can be compared with clinical data of asymptomatic non-obstructive HCM where it has been shown that systolic function in HCM is usually normal or, when abnormal, with high ejection fraction (Wigle et al., 1995).

However, it is worth mentioning that these results are valid in patient suffering from this specific type of HCM only. Indeed, in other cases, HCM may evolve towards ventricular dilatation and wall thinning because of the myocardial fibrosis. The systolic function is then impaired: the ejection fraction is reduced and the end-systolic volume increased.

By analysing the volume variation diagram one can notice that the slope of the volume curves at the filling phase is lower than the one obtained during the simulation of the normal heart. The heart takes more time to relax, the filling is slower. This result is also visible on the video sequence of the simulation where the pathological heart seems to relax slower than the normal organ (see video below). By decreasing the relaxation rate we managed thus to get an anomalous feature consistent with clinical observations (Wigle et al., 1995).

5.3 Conclusion

This experiment constitutes a simulation of hypertrophic cardiomyopathy. Only asymptomatic non-obstructive HCM with global thickening of the left ventricle was considered. Moreover, no focal variations of electro-mechanical features are implemented though very common in patients suffering from HCM. Finally, only two biomechanical parameter, the maximum contraction and the relaxation rate, have been modified according to clinical observation.

Nevertheless, despite these approximations, the 3-D model managed to capture two main features of this type of HCM: a significant increase of the left-ventricular ejection fraction, and a visible impairment of the diastolic function with slower relaxation speed.

For a personnalized simulation, the implementation of focal variations of the electromechanical model of the heart would be needed to better model HCM. In a second stage, valve movements and obstruction may be modelled to further improve the results.

5.4 References

• Maron, B. J., 2002. Hypertrophic Cardiomyopathy. Circulation, vol. 106, pp. 2419-2421
• Wigle, E. D., Rakowski, H., Kimball, B., Williams, W. G., 1995. Hypertrophic Cardiomyopathy, Clinical Spectrum and Treatment. Circulation, No. 92., pp. 1680-1692.

Towards Patient-Specific Heart Models for Decision Support System

1. Current method: From MRI segmentation to 3-D anatomical model of the heart

Tetrahedral mesh of the bi-ventricular myocardium
obtained in (Sermesant et al., 2006)

To adapt the general disease-based model to the patient's anatomy, the most common approach consists in using 3-D imaging to obtain, through image analysis algorithms, a 3-D mesh of the patient's heart.

In (Sermesant et al., 2006) for instance, the 3-D anatomical model used in the disease-based model is created directly from a high-resolution MRI by segmenting the myocardium through classical segmentation techniques (thresholding and mathematical morphology). The marching cube algorithm (Lorensen et al., 1987) is then applied to get a triangulated surface mesh, which is then converted into a volumetric tetrahedral mesh by using the INRIA software GHS3D (see opposite figure).

Similar methods can be applied to the images that are acquired in Health-e-Child. The figures below show an example of volumetric mesh created from a high-resolution MRI of a patient suffering from right-ventricle overload (University College London, Great Ormond Street Children's Hospital) with a supervised algorithm developed at INRIA (Fernandez., 2006). This algorithm allows the segmentation of the four cardiac chambers. The input MRI is first pre-processed through anisotropic filtering to reduce the noise without blurring the edges. Then, geometrical primitives of the cavities are estimated from a preliminary segmentation of the blood pools (via a region-growth method and mathematical morphology). These primitives are then converted into simplex meshes which are automatically deformed according to internal and image-based forces.

Original cardiac MRI acquired at Great Ormond Street Children's Hospital.

Segmentation of the heart from 3-D MRI.
Left panel: Final segmentation. Right panel: A clipping plane is used to show up the myocardium (semi-transparent).
In red: left ventricle, in blue: right ventricle.

Projection of the final mesh into the input MRI. In red, segmentation of the heart.

INTERACTIVE 3D MESH Interactive 3-D mesh. Use the mouse to interact with the mesh

A last method would consist in adapting the bi-ventricular model introduced in the disease-based model to the patient's anatomy. In (Sermesant et al., 2005) for example, the bi-ellipsoidal model of the two ventricles was first generated. Then, the ventricular blood pools were segmented by using boundary- and region-based fuzzy classification (Andriantsimiavona et al., 2003). The resulting segmentation was automatically registered to the model through affine transformation. Finally, local adjustments were performed by using a deformable biomechanical model (Sermesant et al., 2003) and by taking as initial estimation the affine transformation previously computed. The following figures present a result obtained by using this technique.

Bi-ellipsoidal model of heart ventricles (Sermesant et al., 2005).

Adjustment of the model to the patient's anatomy.
(a) Semi-automated segmentation of the ventricular blood pools (white lines).
(b) Affine registration of the segmentation result to the synthetic model.
(c) Local adjustment of the model to the segmentation.

All these methods are nowadays able to extract the heart anatomy from 3-D imaging. The results are satisfying and can be used to simulate the heart function of the patient under study. We plan to make use of such methods to adapt the disease-based model to the patient's anatomy when high-resolution MRI are available.

However, these approaches are sensitive to the quality of the underlying images. If the 3-D image is too noisy or has a slice thickness too high for instance, the segmentation algorithms may have difficulties in extracting the blood pools and the myocardium. These techniques are therefore well adapted to high-resolution MRI or CT but cannot be applied on ultrasonography imaging nor in some 2-D MR images. Consequently, other approaches are required to adapt the bi-ventricular geometrical model to the patient's anatomy when such high-resolution images are not available. Relying on this idea, next section presents preliminary concepts about the generation of anatomical model from sparse 1-D clinical parameters.

References

• Andriantsimiavona, R., Griffin, L., Hill, D., Razavi, R., 2003. Simple cardiac MRI segmentation. International Society for Magnetic Resonance in Medicine Scientific Meeting, vol. 6, p. 951.
• Fernandez, R., 2006. Construction d'un modele anatomique du systeme cardiovasculaire par segmentation d'IRM. Rapport de Stage de Fin d'Etudes.
• Lorensen, W., Cline, H., 1987. Marching cubes: a high resolution 3D surface reconstruction algorithm. Computer Graphics (Proc. SIGGRAPH), vol. 21, no. 4, pp. 163-169.
• Sermesant, M., Forest, C., Pennec, X., Delingette, H., Ayache, N., 2003. Deformable biomechanical models: Application to 4D cardiac image analysis. Med. Image Anal. 7, 475-488.
• Sermesant, M., Rhode, K., Sanchez-Ortiz, G. I., Camara, O., Andriantsimiavona, R., Hegde, S., Rueckert, D., Lambiase, P., Bucknall, C., Rosenthal, E., Delingette, H., Hill, D. L., Ayache, N., Razavi, R., 2005. Simulation of cardiac pathologies using an electromechanical biventricular model and XMR interventional imaging. Med Image Anal, Vol. 9, No. 5., pp. 467-480.
• Sermesant, M., Delingette, H., Ayache, N., 2006. An Electromechanical Model of the Heart for Image Analysis and Simulation. IEEE Trans. Med. Imag., vol. 25, no. 5, pp. 612-625.

Towards Patient-Specific Heart Models for Decision Support System

2. Perspective: From 1-D clinical parameters to 3-D anatomical model of the heart

2.1 Introduction

The second approach that can be used to build a patient-specific heart model is to use 1-D measurements acquired during clinical evaluation instead of 3-D images. This idea comes from the observation that advanced imaging hardware is not systematically available to the clinicians and diagnostic is mostly based on 2-D ultrasonography data. In this way, we propose to integrate into the disease-based heart model 1-D clinical parameters related to the patient's physiology. Simulations would be personalised and hidden parameters that apply to the patient might be extracted without the need of 3-D imaging. Moreover, the generated 3-D patient-specific heart models could be used to improve the accuracy of some measurements, like the volume of the cavities for instance.

As a first step, we rely on geometrical parameters only to define patient-specific heart models. Though mono-dimensional and sparse, these measurements are of great importance to the clinicians since they allow the assessment of the heart anatomy, the efficiency of the myocardium, the cardiac function and, above all, the diagnostic of cardiac pathologies. Moreover, their acquisition is systematic and performed through standardised methods.

Next figure illustrates various measurements that are acquired during the diagnostic process. Size of the ventricles, internal lengths and diameters of the cavity, diameter of the valves and the thickness of the myocardium are the parameters that we plan to use to build the patient-specific heart model.

Clinical measurements used to describe the right ventricle (image from Gaslini).
The white arrows are diameters measured from 2-D or 3-D images.

2.2 Method

The overall algorithm used to generate the patient-specific anatomical model is made up of three main steps:

1. acquisition and input of geometrical parameters,
2. building of the skeleton of the synthetic heart according to these measurements,
3. fitting of deformable surfaces to the skeleton to create the synthetic organ.

To create the 3-D representation of the synthetic organ we make use of the NURBS framework. NURBS, for Non-Uniform Rational B-Spline, is a mathematical model that allows the intuitive generation of sophisticated geometries. Derived from the Bezier's methodology and the B-spline theory, NURBS represents any curves (or surfaces) by:

• a set of control points,
• the related "knots" that specify the way the control points modify the curves (or the surface),
• one (for the curves) or two (for the surfaces) integers that define the order of the NURBS.

Another interesting property of NURBS is the possibility to assign to each control point a weight, which offers more freedom during the generation process.

NURBS surface created using Blender. Left panel: wireframe visualisation of the NURBS surface (in white). The yellow points are the control points with which the surface is modelled. Right panel: solid visualisation of the surface.

Owing to their powerful features, NURBS-based tools can be found in any 3-D modelling software. They are widely used in computer graphics to design sophisticated shapes and are implemented in the most common 3-D modelling suites (3ds Max, Blender...). Advanced techniques based on NURBS have been developed to assist the user. The one which most interests us is the "skinning technique". It consists in modelling the desired shape by defining its skeleton, the software fits automatically NURBS surfaces to this structure to get the final object.

NURBS-based tools are thus suitable for our aims. Blender, a free but powerful 3-D modelling software, have been chosen to generate the patient-specific heart anatomy. Supported by an important community of developers, it implements the essential 3-D modelling tools and offers a Python API for the development of advanced scripts.

A script with a graphical user interface (GUI) has been developed to allow the clinicians to enter the measured parameters and to create automatically the synthetic anatomical heart. In the following the three main steps of the algorithm are detailed.

The first stage of the algorithm consists in providing the clinical measurements. No assumption about the acquisition method has been done and therefore any imaging technique can be utilised, from 2-D ultrasonography to cardiac MRI. The parameters used to generate the synthetic heart have been selected by consensus between the clinicians and the engineers of the Health-e-Child project. The way they are acquired is presented in the diagnostic protocol of the project (Health-e-Child D9.1, 2006) and in the ASE committee recommendations (Lang et al., 2005). Finally, their default values have been set according to the normal values available in the literature (Lang et al., 2005; Ho and Nihoyannopoulos, 2006). The following table lists the parameters required to create the patient-specific heart model along with their default values.

Left ventricle		Right ventricle
Parameters	Default value	Parameters	Default value
Internal length	9.0 cm	Internal length	8.0 cm
Internal diameter	5.0 cm	Internal diameter, major axis	5.0 cm
		Internal diameter, minor axis	3.0 cm
Diameter of the mitral valve	2.8 cm	Diameter of the tricuspid valve, major axis	3.5 cm
		Diameter of the tricuspid valve, minor axis	2.5 cm
Diameter of the aorta valve	1.3 cm	Diameter of the pulmonary valve	1.3 cm
		Diameter of the outflow tract	2.0 cm
		Length of the outflow tract	2.0 cm

Interface used to create patient-specific heart models.

Once the script is provided with the clinical parameters, it creates the skeleton of each ventricle. Because the ventricular geometry cannot be modelled by a single NURBS surface, each ventricle has been divided into two different parts. The first one consists in the "inflow" section of the ventricle, from the atrio-ventricular valve to the apex. The second part corresponds to the "outflow" section, from the arterial valve to the apex. However, this sub-division is transparent to the user, he only has to click on one single button to create the skeleton of a ventricle, Draw LV for the left ventricle and Draw RV for the right ventricle.

Skeleton of the right ventricle generated automatically from clinical parameters

Next, two NURBS surfaces are fitted to the two parts of the skeleton by means of an automated tool that comes up with Blender. Various parameters can be modified to adjust the fitting process but in general their default values offer satisfying results.

NURBS surfaces fitted to the skeleton of the right ventricle

Finally, the two NURBS surfaces are converted into meshes and merged by using an automated Boolean union. Only the external envelop is kept and all the "internal" and artificial walls are removed to get an empty cavity. At last, the resulting mesh can be smoothed to improve the result.

INTERACTIVE 3D MESH

Right ventricle automatically created from clinical parameters.
Left panel: 2-D view of the 3-D mesh.
Right panel: interactive 3-D mesh (use the mouse to interact with the 3-D model)

2.3. Conclusion

Relying on the observation that 3-D imaging is not systematically available, a preliminary algorithm that automatically creates 3-D heart geometries from well-known and standardised 1-D clinical parameters has been presented. The underlying principle are promising and clinicians involved in the Health-e-Child project have shown a great interest in the method.

The exploratory algorithm developed so far relies on geometrical parameters only and on the NURBS framework to model the cavities. Blender software has been used for the implementation and a graphical user interface has been developed.

Though NURBS-based tools for the modelling of the heart anatomy are available in the literature, the overall majority, to our knowledge, depends on the manual definition of landmarks in the image to build the NURBS surfaces and generate the synthetic heart model (Segars et al., 1999). Our approach is thus distinct from these methods since no landmarks are used and the shape is controlled by clinical and meaningful measurements.

2.4. References

• Health-e-Child Integrated Project, 2006. Report on diagnostic coding system and data collection protocols.
• Ho, S. Y., Nihoyannopoulos, P., 2006. Anatomy, echocardiography, and normal right ventricular dimensions. Heart, Vol. 92, pp. 2-13.
• Lang, R. M., Bierig, M., Devereux, R. B., Flachskampf, F., A., Foster, E., Pellikka, P., A., Picard, M., H., Roman, M., J., Seward, J., Shanewise, J., S., Solomon, S., D., Spencer, K., T., Sutton, M., St J., Stewart, W., J., 2005, Recommendations for Chamber Quantification: A Report fromt the American Society of Echovardiography's Guidelings end Standards Committee and the Chamber Quantification Writing Group, Developed in Conjunction with the European Association of Echocardiography, a Brandch of the European Society of Cardiology. J. Am. Soc. Echocardiogr, no 18, pp. 1440-1463
• Segars, W. P., Lalush, D. S., Tsui, B. M. W, 1999. A realistic spline-based dynamic heart phantom. Nuclear Science, IEEE Transactions on, Vol. 46, No. 3., pp. 503-506.

Personalised Simulation of Myocardium Electromechanics
and Pulmonary Valve Replacement Surgery in
Repaired Tetralogy of Fallot: a Case Study

1. Introduction: Medical Background and Challenges

Tetralogy of Fallot is a congenital heart disease that mainly affects the right ventricle (denoted as RV in the following). It occurs in about 1 child out of 2500 (Hoffman et al., 2002) and it is characterised by four heart defects:

• Inter-ventricular septal defect
• Pulmonic stenosis
• Hypertrophic right ventricle
• Overriding aorta

Comparison between a healthy heart and Tetralogy of Fallot (from Wikipedia)

Because of these defects, Tetralogy of Fallot requires surgical repair early in infancy. Mortality and morbidity of this initial intervention is nowadays very low but the repair may lead to the destruction of the pulmonary valves. Regurgitations can appear, resulting in extreme dilation of the right ventricle and loss of cardiac function later during the child life. Pulmonary valve may thus be replaced to retrieve a normal cardiac function.

Determining the best time for re-intervention is a clinical challenge. Because the child is growing, replacing the valves too early may not be necessary and may result to be inefficient on the long term. Besides, valves lifetime is about 10 years. Hence, the earlier they are replaced, the earlier they must be changed, resulting in repeated interventions. On the other hand, waiting too much may result fatal. The extreme dilation of the right ventricle may cause arrhythmia or other electro-physiological troubles. Furthermore, if too dilated, the right ventricle may not be able to recover a normal function even after valve replacement.

Several pulmonary valve replacement strategies are available. In Health-e-Child, we are interested in two competing approaches:

1. Pulmonary valve replacement alone. Usually, this pulmonary valves are replaced surgically but minimally-invasive techniques have been proposed recently (Khambadkone et al., Circulation 2005). With this approach, the heart adapts itself to its new conditions owing to natural remodelling.
2. Pulmonary valve replacement and direct right ventricle remodelling. The surgeon not only replaces the valves but also resects any scars and fibrosis regions of the myocardium to improve the cardiac function and reduce the right ventricle volume (del Nido et al., 2006).

Nevertheless, choosing which therapy to apply for a given patient is nowadays complex. The effects of these therapies upon the cardiac function are difficult to predict. No clinical evidences nor rules have been found to help the cardiologists in taking the right decision.

Based on these clinical rationales, we investigate in Health-e-Child how cutting-edge image analysis algorithms and personalised electromechanical models of the heart can help in evaluating the cardiac function of a patient and can predict the effects of right pulmonary valve replacement therapies. To that purpose, we propose a three-step approach:

• Step 1: Building a patient-specific anatomical model of the heart from clinical gated cine-MRI.
• Step 2: Simulating patient cardiac function by personalising an electromechanical model of the heart
• Step 3: Simulating the two PVR therapies using soft-tissue intervention platform

The pipeline is tested on a young Health-e-Child patient (age=17). Details and first results of the experiments are presented in the following sections.

References

• Khambadkone, S., Coats, L., Taylor, A., Boudjemline, Y., Derrick, G., Tsang, V., Cooper, J., Muthurangu, V., Hegde, S. R., Razavi, R. S., Pellerin, D., Deanfield, J., and Bonhoeffer, P.: Percutaneous Pulmonary Valve Implantation in Humans: Results in 59 Consecutive Patients, Circulation 112(8), 1189–1197, 2005
• Hoffman, J.I.E., and Kaplan, S.: The incidence of congenital heart disease, Journal of the American College of Cardiology 39(12), Elsevier, 1890-1900, 2002
• del Nido, P. J.: Surgical Management of Right Ventricular Dysfunction Late After Repair of Tetralogy of Fallot: Right Ventricular Remodeling Surgery, Semin Thorac Cardiovasc Surg Pediatr Card Surg Annu. 9(1), 29-34, 2006

Personalised Simulation of Myocardium Electromechanics
and Pulmonary Valve Replacement Surgery in
Repaired Tetralogy of Fallot: a Case Study

2. Creating a Personalised Anatomical Model of the Heart

Simulating patient cardiac function cannot be performed without considering the geometry of patient heart. Therefore, we need to extract from clinical gated cine-MRI a reliable representation of the biventricular myocardium of the patient. Furthermore, having 3D+t delineation enables us estimating dynamic measurements such as blood pool volume variations, which help us in personalising the electromechanical model.

The anatomical model comprises several information regarding the patient heart. First, it represents the shape of the organ. This can be achieved by segmenting the myocardium from clinical cine-MRI. Next, it defines the anatomical regions that are characterised by different motion patterns, for regional personalisation of the electromechanical model. Finally, it encodes the direction of the myocardium fibres to simulate the anisotropic properties of the myocardium.

2.1. Delineating Myocardium Shape

Numerous methods have been proposed in the literature to efficiently segment the heart from cine MRI (see Part III). However, these methods are not readily adapted to our application. In children, myocardium is thinner, images are less contrasted and motion artifacts are more frequent. Moreover, repaired tetralogy of Fallot leads to abnormal myocardium motion patterns, extreme dilation and tremendous variability of the right ventricle shape. Using atlas-based methods is thus not appropriate, these approaches may fail because of the paediatric image specificities. To tackle these challenges, we devised specific tools that requires few manual interactions.

Delineation of the right-ventricle endocardium

In Health-e-Child, the right-ventricle endocardium is segmented on all the frames of the cardiac sequence by fitting an anatomically accurate geometrical model. Its position, orientation and scale in the images are determined using minimal user interaction. Boundaries are locally adjusted by training a probabilistic boosting tree classifier with steerable features (Zheng et al., ICCV 2007). The resulting delineation is then tracked throughout the cardiac cycle using an optical flow method.

Delineation of the left-ventricle endocardium and of the epicardium

Interactive delineation of the left endocardium. The surface (in wireframe) is modelled in real-time by interactively placing control points inside (green point), on (red points) and outside (blue points) the surface.

The left-ventricle endocardium and the epicardium are delineated on the first frame of the gated cine MRI by using an interactive tool based on variational implicit functions (Turk et al., Technical Report 1999). The user places control points inside, on and outside the desired surface. The algorithm computes in real-time the implicit function that interpolates those points and extracts its 0-level set (Toussaint et al., VCBM 2008).

Once the desired surface is modelled, diffeomorphic non-linear registration (Vercauteren et al., MICCAI 2007) is applied to automatically propagate the surfaces throughout the cardiac cycle.

Construction of the Dynamic Bi-Ventricular Myocardium

The surfaces computed in the previous steps are combined together to get the dynamic segmentation of the myocardium. To that aim, we combine the binary masks of the surfaces, thus obtaining a dynamic mask of the biventricular myocardium. The valve plane is manually defined using CardioViz3D and muscle consistency is ensured by preserving a minimal thickness of 3 mm (mean thickness of a healthy right ventricle myocardium).

Next, we extract the 3D surface of the myocardium on the first frame. Finally, simplex-based deformable surfaces (Montagnat et al., MedIA 2004) are used to propagate the myocardium mesh throughout the cardiac cycle.

Quantitative Measurements

The dynamic contours we have obtained enable us to compute quantitative and reproducible measurements related to the cardiac function. In particular, blood pool volumes can be measured at each time frame. For the patient we are considering, ejection fraction is found normal for the left ventricle (61%) but slightly below the normal values for the right ventricle (41%). However, right ventricle volume is twice as great as left ventricle volume (see curve below).

Additionally, the 3D representation of the beating heart enables us to estimate 3D indexes that are not readily available in clinical routine. For instance, radial displacements of the myocardium, computed at each vertex of the 3D mesh, highlight the abnormal motion of the right-ventricle outflow tract. This region dilates when the healthy myocardium contracts. Such observation is of crucial interest when personalising the electromechanical model. We can indeed locally personalise this region by choosing appropriate parameters in order to recover the pathological motion.

Volume variations throughout the cardiac cycle of the left (in red) and right (in blue) ventricles

Get the Flash Player to see this player. Radial displacement of the myocardium. Observe the dyskinetic motion of the right ventricle outflow tract due to the pathology

2.2. Personalised Myocardium Anatomy

The anatomical model can now be defined. From the volume curves, we select the mid-diastole time frame. The related surface mesh is transformed into a tetrahedral volume mesh. From the endocardia and epicardium surfaces, left and right ventricles are automatically detected and defined in the volumetric mesh. At last, the pathological region is specified in the anatomical model to simulate its pathological motion.

2.3. Myocardium Fibres

Myocardium fibres play a crucial role in cardiac biomechanics. Unfortunately, measuring fibre direction in-vivo is still difficult. To tackle this limitation, we use a computational model based on observations done on anatomical dissections or post-mortem diffusion tensor images. These studies showed that fibre orientation varies from -70° on the epicardium to 0° at mid-wall to +70° on the endocardium (Arts et al., 2001, Peyrat et al., TMI 2007, Sanchez-Quintana et al., 1996). Synthetic fibres are thus created by linearly interpolating their orientation with respect to the short axis plane, from -90° on the epicardium to 0° at mid-wall to +90° on the endocardium. The angles are slightly over-estimated to account for rasterisation effects.

Anatomical model of the patient heart. For regional personalisation of the electromechanical model, we define in red the left ventricle, in orange the right ventricle and in white the observed dyskinetic area.

Myocardium fibres estimated using a computational model adapted to the pathology (Sanchez-Quintana et al., 1996)

2.4. Conclusions

From the cine-MRI of the patient, we have semi-automatically delineated the biventricular myocardium. Then, from this segmentation, we have built an anatomical model of the patient heart, which represents the shape of the organ, defines some regions of interest in view of regional personalisation of the electromechanical model, and sets the direction of the myocardium fibres. Furthermore, the dynamic segmentation has provided us with volume curves and radial displacements, measurements which will help us in personalising the model, as described in the next section.

References

• Arts, T., Costa, K. D., Covell, J. W., and McCulloch, A. D.: Relating myocardial laminar architecture to shear strain and muscle fiber orientation, Am J Physiol Heart Circ Physiol 280(5), 2001
• Montagnat, J., and Delingette, H.: 4D Deformable Models with temporal constraints: application to 4D cardiac image segmentation, Medical Image Analysis 9(1), 87-100, 2005
• Peyrat, JM., Sermesant, M., Pennec, X., Delingette, H., Xu, C., McVeigh, E. R., and Ayache, N.: A Computational Framework for the Statistical Analysis of Cardiac Diffusion Tensors: Application to a Small Database of Canine Hearts, IEEE TMI 26(11), 1500-1514, 2007
• Sanchez-Quintana, D., Anderson, R. H., and Ho, S. Y.: Ventricular myoarchitecture in tetralogy of Fallot., Heart 76(3), 280-286, 1996
• Toussaint, N., Mansi, T., Delingette, H., Ayache, N., and Sermesant, M.: An Integrated Platform for Dynamic Cardiac Simulation and Image Processing: Application to Personalised Tetralogy of Fallot Simulation, Proc. Eurographics Workshop on Visual Computing for Biomedicine (VCBM), 2008
• Turk, G., and O'Brien, J.: Variational Implicit Surfaces, Georgia Institute of Technology, 1999
• Vercauteren, T., Pennec, X., Perchant, A., and Ayache, N.: Non-parametric Diffeomorphic Image Registration with the Demons Algorithm, Proc. MICCAI 2007
• Zheng, Y., Barbu, A., Georgescu, B., Scheuering, M., and Comaniciu, D.: Fast Automatic Heart Chamber Segmentation from 3D CT Data Using Marginal Space Learning and Steerable Features, Proc. ICCV 2007, 1-8, 2007

Personalised Simulation of Myocardium Electromechanics
and Pulmonary Valve Replacement Surgery in
Repaired Tetralogy of Fallot: a Case Study

3. Personalising an Electromechanical Model to Simulate Patient Cardiac Function

Once the personalised anatomy is built, we can personalise the electromechanical model presented in Part II, section 1, to simulate patient cardiac function. This model is made up of three main building blocks: electrophysiology, to simulate the electrical command; biomechanics, to model the pump function; and boundary conditions, to place the virtual heart in the body. We must personalise each one of these elements to simulate patient cardiac function.

3.1. Electrophysiology

Model

Simulating electrophysiology is important since tetralogy of Fallot patients may suffer from electrophysiological troubles, which can lead to sudden death. In particular, bundle-branch blocks are frequent, mainly due to the initial surgical repair and subsequent fibrosis. QRS duration can also be significantly longer than in healthy subjects.

Because we aim at simulating these abnormalities but also at personalising the model using clinical data, we require a model that is accurate but simple enough to be controlled by few clinical parameters. In Health-e-Child, we use an anisotropic multi-front Eikonal model (Sermesant et al., FIMH 2007). This model provides, in non-extreme cases, similar results as those obtained with more complex models. However, it is controlled by only two clinically meaningful parameters and it is significantly faster (of the order of the seconds). The parameters to be set are:

• Electrical conductivity, which can be estimated using electrocardiograms or catheter measurements,
• Electrical anisotropy. Nominal values can be found in experimental studies.

Personalisation

For our patient, only the electrocardiogram is available, where no evident abnormality is observed. Nominal values available in the literature are thus used: conductivity is set to 500 m/s and anisotropy to 3. As illustrated by the video on the left, the electrical wave is initialised at the end of the Purkinje bundles, near the apex.

3.2. Biomechanics

Model

Once the electrical command personalised, we must personalise the cardiac contraction to simulate patient heart motion. As presented in Part II, section 1, the biomechanical model is constituted of two main elements: a geometrical model and a biomechanical constitutive law that models the properties of the myocardium tissue (Sermesant et al., TMI 2006).

Regarding the geometrical model, we rely on a tetrahedral representation of the organ. This enables us to use the personalised anatomies described in the previous section. Linear elasticity is used to model cardiac motion.

The constitutive law used in Health-e-Child is based on a multi-scale model, from nano-motors to the organ (see Part II, section 1). This model is controlled by the following parameters:

• Tissue mass: nominal values can be found in the literature and can be reasonably assumed constant.
• Anisotropy: nominal values can be found in the literature and can be reasonably assumed constant.
• Elasticity: nominal values available in the literature but can be qualitatively estimated from cine-MRI or elastography.
• Contractility: nominal values available in the literature but can be qualitatively estimated from cine-MRI or ultrasound images.

Personalisation

The biomechanical model is personalised by alternatively modifying the elasticity and contractility parameters based on observations done from the cine-MRI and the dynamic segmentations. The other parameters are set according to the values found in the literature (mass=1070kg/m3, anisotropy=3). Elasticity and contractility are set independently for three myocardium regions: left ventricle, right ventricle and dyskinetic area. A trial-and-error strategy is carried out to estimate these parameters. Cardiac function is simulated and compared to the cine MRI, 3D segmentation and other clinical observations such as 2D-US end-systolic pressures. Then, contractility is increased or decreased, to retrieve the measured ejection fractions and blood pool volume variations. Elasticity is modified to calibrate the simulated motion to the observed motion. Finally, the contractility of the dyskinetic area is set to 0 to simulate the pathology.

3.3. Boundary Conditions

Model

The last step consists in placing the virtual heart in his environment, i.e. the body. We simulate the four cardiac phases: contraction, ejection, relaxation and filling (Sermesant et al., TMI 2006). Atrial and arterial pressures are also taken into considerations. Finally, pulmonary regurgitations are simulated using a simple model that relies on data measured from 2D ultrasound. The measured regurgitation flows are directly integrated in the virtual heart (Mansi et al., 3DPH 2008).

Personalisation

Regurgitation flows are measured using 2D ultrasound Doppler images and directly integrated into the model. Doppler measurements also provide us with the end-systolic right-ventricule pressure (15mmHg), which we use to personalise the pulmonary artery pressure. Finally, the atrial pressures are set to 5mmHg for the right atrium and 12mmHg for the left atrium. End-systolic aorta pressure is set to 75mmHg following Health-e-Child cardiologists suggestions.

3.4. Results

Several quantitative indexes are obtained from the personalised simulation. First, simulated blood pool volumes are computed. Of course, this data was already available from the segmentations but simulated volumes enable us to calibrate the model. In that way, the biomechanical parameters (elasticity and contractility) are set such that the simulation agrees with the observed volumes and ejection fractions.

However, the key advantage of performing personalised simulations is that parameters that are not readily available in clinical routine can be simulated. For instance, we can simulate the variation of the ventricular, atrial and arterial pressures during the cardiac cycle. We are also able to draw pressure-volume loops, which provide precious information about myocardium efficiency.

Simulated volume variations (right ventricle in blue, left ventricle in red) and segmented volume variations (in green). The electromechanical model of the heart managed to recover the global cardiac function of the patient. Vertical bars limit the four cardiac phases.	Simulated pressures in the left ventricle (solid red curve), in the right ventricle (solid blue curve), in the aorta (dashed red curve) and in the pulmonary artery (dashed blue curve)
Pressure-volume diagram deduced from the two previous diagrams. Observe the effect of the pulmonary regurgitations during the "iso-volumetric" phases

Simulated contractions and 3D strain can also be assessed, to understand what are the effects upon the cardiac function of local abnormalities such as the dyskinetic region we observed in the images. Finally, radial displacements can be calculated and compared to the ground truth to refine the personalisation.

3.5. Conclusions

We now have a personalised electromechanical model of the heart that qualitatively simulates the cardiac function of our patient. The virtual heart can now be used to simulate parameters that are difficult to obtain in clinics. But more importantly, the parameters can be modified to simulate the effects of a therapy, as we will describe in the following section.

References

• Mansi, T., Andre, B., Lynch, M., Sermesant, M., Delingette, H., Boudjemline, Y., and Ayache, N.: Virtual Pulmonary Valve Replacement Interventions with a Personalised Cardiac Electromechanical Model, Workshop on 3D Physiological Human (in press), December 2008
• Sermesant, M., Delingette, H., and Ayache, N.: An Electromechanical Model of the Heart for Image Analysis and Simulation, IEEE TMI 25(5), 612-625, 2006
• Sermesant, M., Konukoglu, E., Delingette, H., Coudiere, Y., Chinchapatnam, P., Rhode, K., Razavi, R., and Ayache, N.: An Anisotropic Multi-front Fast Marching Method for Real-Time Simulation of Cardiac Electrophysiology, Proc. FIMH 2007, Proc. FIMH 2007, 160-169, 2007

Personalised Simulation of Myocardium Electromechanics
and Pulmonary Valve Replacement Surgery in
Repaired Tetralogy of Fallot: a Case Study

4. Simulating Pulmonary Valve Replacement Therapies

4.1. Medical Context

In the previous sections, we have described how we personalise an electromechanical model of the heart to simulate patient cardiac function. We now have a virtual heart that represents the pathophysiology of the patient. We can use this tool to simulate therapies, in particular we can test the two pulmonary valve replacement (PVR) strategies presented in the introductory section.

Pulmonary valve replacement is becoming a key therapy for patients suffering from chronic pulmonary valve regurgitations, and in particular for patients suffering from repaired Tetralogy of Fallot. The therapy consists in replacing new valves between the right-ventricle outflow tract and the pulmonary artery to remove the regurgitations. Nowadays, two competing approaches are available:

1. Valve replacement alone. Only the valves are replaced. If the right ventricle is too dilated, the surgeon may have to reduce the diameter of the right-ventricle outflow tract to place the valves, but no extra surgery is performed. The heart adapts itself to its new condition by natural remodelling.
2. Valve replacement with direct right-ventricle remodelling. The surgeon replaces the valves and he directly resects any myocardium regions that are impaired by the pathology. Scars and fibrosis are removed to improve the efficiency of the cardiac muscle. Right ventricle is explicitly reduced (del Nido et al., 2006).

However, the effects of these therapies on the cardiac function of the patient cannot be predicted easily. No clinical evidences are available to help the cardiologist in choosing the appropriate therapy. Therefore, we investigate in Health-e-Child how the virtual heart can be used to simulate the PVR therapies and to predict their effects on the cardiac function. The therapies are simulated as follows:

1. Valve replacement. Because replacing the valves only affects the regurgitations, we simulate this intervention by only disabling the regurgitations in the models.
2. Direct remodelling of the right-ventricle. Virtual soft-tissue intervention is carried out using SOFA framework as explained below.

4.2. Virtual Surgery

In our example, virtual surgery consists in removing the dyskinetic area we have observed in the images. As suggested by the biomechanical parameters we found when personalising the model, this area presents with a very low contractility, probably due to fibrosis or to the aneurysm. By removing this impaired area we hope to improve the efficiency of the myocardium.

More precisely, the virtual surgery is made up of three interactive and real-time steps:

1. Resecting the impaired myocardium tissue: we remove the aneurysm and the surrounding tissue that might be affected by fibrosis. This task is carried out through a point-and-click approach. The user first points at the center of the area to resect, then at the boundary, and all the mesh tetrahedra that lay within a sphere defined by the two points are removed (see movie below).
2. Remodelling the right ventricle: the free-wall is smoothly and carefully brought near the septum, to close the space resulting from the resections. The user clicks to grab the muscle and moves it. Simulation time-step can be modified in real-time to enable smooth and controlled displacements. Indeed, this operation is very sensitive since the virtual heart is empty of blood, which makes the remodelling complex. Furthermore, the final shape of the organ must be anatomically plausible despite of the resections.
3. Suturing the myocardium tissue: we sew the myocardium to definitively close the right-ventricle. After suture, the scar is defined to simulate the possible effects of the surgical intervention in that region, where electrical conductivity and myocardium contractility may result damaged by the intervention.

Virtual RV volume reduction surgery. From left to right: original mesh, after resection, during attachment, final mesh. Colour lines: fibre orientations. Black area: postoperative scar.

4.3. Results

Once the intervention is performed on the anatomical model, we run the electromechanical simulation using the same parameters as those estimated during the personalisation, except for the regurgitations, which are disabled (to simulate the replacement of the valves).

Volumes and pressures are computed at each time step. From the results, we can conclude that, for this patient, the best therapy would be replacing the valves and directly remodelling the myocardium. With this strategy indeed, the right ventricle volume is reduced, as a direct consequence of the intervention, but the right ejection fraction is significantly higher. Furthermore, left ventricle function is also improved, which highlights the relationship between the two ventricles. Finally, the right ventricle pressure at end-systole slightly decreases because of the scar, which modifies the myocardium motion in this area and decreases its efficiency.

RV volume curves and pressure-volume loops: (solid green) from segmentation; (solid black) preoperative simulation; (dashed blue) simulation after PPVR; (dash-dotted red) simulation after PVR with volume reduction. Vertical bars delineate the cardiac phases. Regurgitations are visible during the isovolumetric phases, when volume should be constant. As expected, RV volume reduction resulted in a shift of the curves.

4.4. Conclusions

For this patient, we found that valve replacement with volume reduction would provide a better immediate postoperative outcome. Probable reasons may be the removal of the dyskinetic area and the direct reduction of the RV volume. However, this procedure is very invasive and may be hazardous for the patient, with possible postoperative side-effects such as electro-physiological troubles due to the surgical scar. On the other hand, the effects of valve replacement alone are often visible late after replacement, after the myocardium adapts itself to its new loading conditions. A model of the natural remodelling of the myocardium is thus important to simulate the long-term postoperative effects of valve replacement.

References

• del Nido, P. J.: Surgical Management of Right Ventricular Dysfunction Late After Repair of Tetralogy of Fallot: Right Ventricular Remodeling Surgery, Semin Thorac Cardiovasc Surg Pediatr Card Surg Annu. 9(1), 29-34, 2006

Personalised Simulation of Myocardium Electromechanics
and Pulmonary Valve Replacement Surgery in
Repaired Tetralogy of Fallot: a Case Study

5. Conclusions

We have studied the potential of image processing techniques, electromechanical models and virtual soft-tissue intervention platforms to perform virtual and personalised assessment of pulmonary valve replacement (PVR) therapies. The results have been promising and suggested that such tools might be used by the clinicians to test different PVR therapies.

A modular framework has been proposed. The electromechanical model presented in the previous sections has been adjusted to simulate the cardiac function of an enrolled patient who suffers from repaired Tetralogy of Fallot. An interactive soft-tissue intervention platform has been adapted to simulate right-ventricle remodelling. The framework managed to realistically simulate the cardiac function of the selected patient. For this patient, we found that PVR with RV volume reduction would yield significative improvements just after the intervention.

Owing to the modularity of our framework, more elaborate models can now be used. In particular, future works include: improvements of the regurgitation model, to take into account the other cardiovascular compartments, modelling of the long-term effects of the PVR therapies, by considering myocardium remodelling, automated estimation of biomechanical parameters, comprehensive validation using preoperative and post-operative clinical data. Such in-silico surgery simulation opens up possibilities in therapy planning, for a more predictive and personalised medicine.

iLogDemons, a Consistent Incompressible Image Registration Algorithm

1. Introduction: Medical Background and Challenges

Tissue tracking in sequences of medical images is an important task in many applications, either for therapy guidance or diagnosis. However there is no easy way to achieve it, even interactively. A standard approach is now to use non-linear image registration to estimate the spatial transformation between different images, for instance brain shift (Ferrant et al., 2001). However, for such clinical applications one-to-one mapping must often be ensured and the model of transformations that is used must be consistent with the properties of the tissue to track. In particular, constraining the registration to be volume-preserving showed great improvement when tracking incompressible tissues such as the heart (Bistoquet et al., 2008; Mansi et al., 2009).

Mathematical frameworks based on diffeomorphic deformations (Beg et al., 2005; Bossa et al., 2007) have been developed to get one-to-one mappings between the images to register. Among them, logDemons (Vercauteren et al., 2008) is an efficient demons-based registration algorithm that provides diffeomorphic transformations parameterised by stationary velocity fields. However, although mathematical justifications of demons optimisation have been provided (Vercauteren et al., 2008), theoretical foundations of the Gaussian regularisation still has to be consolidated (Pennec et al., 1999), which makes the algorithm difficult to adapt to other deformation models.

Several approaches for incompressible image registration have been proposed. A first method is to constrain the Jacobian determinant of the transformations to equal one. However, this constraint is computationally demanding due to its non-linearity (Rohlfing et al., 2003). Linear approximations have been proposed (Bistoquet et al., 2008) but volume drifts may appear at large deformations. Velocity fields can be made incompressible by constraining them to be divergence-free. Thus, incompressible fluid registration is achieved by projecting the update velocity onto the space of divergence-free vector fields (Saddi et al., 2007). Nevertheless, the fluid model might not be appropriate for tracking biological soft tissues like myocardium.

In Health-e-Child we developed an efficient and consistent framework for demons-based incompressible registration. We first propose a mathematical justification of the Gaussian regularisation, which enables to integrate the incompressibility constraint seamlessly by working on the space of divergence-free velocity fields. The main advantages of this are:

• the constraint is linear with little computational overhead
• the parameter of the deformations are constrained: no volume drifts appear
• the transformation minimises a constrained energy functional: the optimal incompressibility field according to the logDemons minimisation scheme is found. The algorithm was validated against synthetic data and applied on clinical cardiac MRI to estimate 3D myocardium motion.

1. Introduction: Medical Background and Challenges

Background: Log-Domain Diffeomorphic Demons

LogDemons algorithm estimates a dense non-linear transformation φ that best aligns a template image T to a reference image R (Vercauteren et al., 2008). φ belongs to the space generated by the one-parameter subgroups of diffeomorphisms. They are parameterised by stationary velocity fields v through the exponential map φ=exp(v) (Arsigny et al., 2006). The images R and T are registered by minimising in the space of velocities, called log-domain, the energy functional:

where σ_i² relates to the noise in the images and σ_d² controls the regularisation strength. This equation is function of two variables: The velocity field v parameterises the transformation φ to recover whereas v_c parameterises an intermediate transformation φ_c = exp(v_c) that models the correspondences between the voxels of the two images. During the optimisation step, E(v,v_c) is minimised with respect to v_c. This amounts to finding the optimal matching between R and T without considering the regularisation. Under the diffeomorphic update rule φ_c ← φ o exp(Δv), the optimal update velocity writes

J(x) is the symmetric gradient J(x) = (∇R(x) + ∇(T o φ)(x)) / 2. The correspondence velocity v_c is then updated using the first order approximation of Baker-Campbell-Hausdorff (BCH) formula:

where the Lie bracket [·,·] is defined by [v,Δv] = (∇v) Δv - (∇Δv)v.

This approach has experimentally shown promising results in terms of image registration and statistics on diffeomorphisms (Bossa et al., 2007; Vercauteren et al., 2008). Finally, the regularisation step estimates the optimal regularised transformation φ by minimising E(v,v_c) with respect to v, which is approximated by smoothing the correspondence velocity v_c with a Gaussian kernel G_σ. However, how G_σ relates to E(v,v_c) remains to be consolidated (Pennec et al., 1999).

Revisiting Demons Gaussian Regularisation

A consistent mathematical formulation of demons regularisation is required to integrate incompressibility in logDemons. In (Nielsen et al., 1994), the authors demonstrate that Gaussian filtering solves the Tikhonov estimation problem with equal weighting of the spatial derivatives in the Taylor series sense. We thus replace the logDemons regulariser ||∇v||² by the Tikhonov regulariser to get:

In this equation, Ω is the image domain and ∂_ik..il denotes the composition of spatial derivatives ∂_ik..∂_il. The parameter σ_x² has been introduced into the regulariser to simplify calculations. More importantly, the regularisation weight σ_d² is now function of the derivative orders to preserve the shape of the impulse response related to the regulariser (Nielsen et al., 1994). The previous equation is minimised by linearising its first term using the zeroth order approximation of BCH formula, log(φ^-1 o φ_c) = v_c - v. The resulting equation, which is exactly a Tikhonov estimation problem, is minimised in the Fourier domain. The optimal velocity field v verifies

which is exactly demons Gaussian smoothing G_σ, with σ² = 2/σ_d². The width of the Gaussian kernel corresponds to the strength of the regularisation.

This formulation is built up on two key elements. First, the correspondence field φ_c decouples the regularisation from the optimisation, making the minimisation of the registration energy independent of the optimisation. Second, the coupling term ||log(φ^-1 o φ_c)||² is approximated as a least-square problem to get the Gaussian smoothing. These two elements must be ensured in any demons algorithm to justify the Gaussian regularisation.

Incompressible LogDemons

The proposed regulariser enables to integrate the incompressibility constraint into the algorithm. A transformation φ is locally incompressible if its Jacobian determinant det(∇φ)=1. In fluid dynamics, one uses the infinitesimal version: a fluid is said incompressible if the divergence of its velocity v is null. For diffeomorphic transformations one can show that the converse is true: integrating divergence-free velocities over time yields incompressible deformations (Evans, 1998). Making LogDemons incompressible thus consists in constraining the velocity fields v to be divergence-free. Demons optimisation step is not modified, as it optimises v_c only, but demons regularisation energy is now optimised under the divergence-free constraint, which amounts to minimising the Lagrange function:

In this equation, the Lagrange multiplier p is a scalar function of the Sobolev space H₀¹(Ω) that vanishes at infinity. Optima of the previous equation are found by solving ∂_vP(v,p)=0:

with p=0 at the domain boundaries ∂Ω. The divergence of the previous equation, under the optimal condition ∇·v=0, yields Δp = ∇· v_c with 0-Dirichlet boundary conditions, which can be solved independently of v to get p. g = v_c - ∇p is thus the L₂ projection of v_c onto the space of divergence-free vector fields. From the previous section we deduce the optimal incompressible velocity field: v ←G_σ⊗g.

Particular care must be taken when incompressibility is required within a subdomain Γ ⊂ Ω only, for tracking localised incompressible organs like the heart. This is achieved by defining p ∈ H₀¹(Γ), p=0 on Ω/Γ. Although Gaussian smoothing theoretically preserves the divergence, in practice unconstrained velocities close to Γ boundaries (∂Γ) may leak inside the incompressible domain due to the Gaussian convolution, ultimately resulting in volume drifts. Yet, vector derivatives and well-designed Gaussian filters, like Deriche recursive filters, commute. We therefore replace the theoretical "project-and-smooth" strategy by a "smooth-and-project" approach that preserves divergence. To limit numerical instabilities, a smooth domain transition is implemented in a narrow band around Γ by diffusing the pressure field p using heat-transfer equation (Evans, 1998). The following algorithm summarises the main steps of our method, henceforth termed iLogDemons, which was implemented using ITK (Dru et al., 2009) and PETSc.

References

1. Arsigny, V., Commowick, O., Pennec, X., Ayache, N.: A log-euclidean framework for statistics on diffeomorphisms. In MICCAI'06, LNCS, v. 4190, pp. 924 Springer (2006)
2. Beg, M., Miller, M., Trouvé, A., Younes, L.: Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms. IJCV 61(2), 139-157 Springer (2005)
3. Bistoquet, A., Oshinski, J., Skrinjar, O.: Myocardial deformation recovery from cine MRI using a nearly incompressible biventricular model. Medical Image Analysis 12(1), 69-85 Elsevier (2008)
4. Bossa, M., Hernandez, M., Olmos, S.: Contributions to 3D diffeomorphic atlas estimation: application to brain images. In MICCAI'07, LNCS, vol. 10, p. 667 Springer (2007)
5. Dru, F., Vercauteren, T.: An ITK implementation of the symmetric log-domain diffeomorphic demons algorithm. Insight Journal (May 2009)
6. Evans, L.C.: Partial Differential Equations. (June 1998)
7. Ferrant, M., Nabavi, A., Macq, B., Jolesz, F., Kikinis, R., Warfield, S.: Registration of 3D intraoperative MR images of the brain using a finite element biomechanical model. IEEE TMI 20(12), 1384-1397 (2001)
8. Hinkle, J., Fletcher, P., Wang, B., Salter, B., Joshi, S.: 4D map image reconstruction incorporating organ motion. In IPMI'07, LNCS, p. 687. Springer (2009)
9. Mansi, T., Peyrat, J.M., Sermesant, M., Delingette, H., Blanc, J., Boudjemline, Y., Ayache, N.: Physically-constrained diffeomorphic demons for the estimation of 3D myocardium strain from cine-MRI. In FIMH'09, LNCS, v. 5528, p. 201-210 Springer (2009)
10. Nielsen, M., Florack, L., Deriche, R.: Regularization, scale-space, and edge detection filters. In JMIV, 291-307 Springer (1997)
11. Pennec, X., Cachier, P., Ayache, N.: Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent. In MICCAI'99. LNCS, vol. 1679, pp. 597-605. Springer (1999)
12. Rohlfing, T., Maurer Jr, C., Bluemke, D., Jacobs, M.: Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint. IEEE TMI 22(6), 730-741 (2003)
13. Saddi, K.A., Chefd'hotel, C., Cheriet, F.: Large deformation registration of contrast-enhanced images with volume-preserving constraint. In: SPIE Medical Imaging. vol. 6512. (2007)
14. Vercauteren, T., Pennec, X., Perchant, A., Ayache, N.: Symmetric log-domain diffeomorphic registration: A demons-based approach. In MICCAI'08, LNCS, vol. 5241, pp. 754-761. Springer (2008)
15. Moore,C.,Lugo-Olivieri,C.,McVeigh,E.,Zerhouni,E.: Three-dimensional systolic strain patterns in the normal human left ventricle: Characterization with tagged MR imaging. Radiology 214(2), 453-466 (2000)

iLogDemons, Experiments on Synthetic Data

iLogDemons were tested on synthetic datasets with known ground truth. Eight 3D volume-preserving whirl transformations φ_α were created with whirl angles α= 10° to 80° (Saddi et al., 2007). Within the whirl domain, the L₂-norm varied from 0.52mm to 4.78mm while the Jacobian determinant stayed close to one (worse value |∇φ_α=80°|= 1±0.04 (mean ± standard deviation SD)). A 3D isotropic Steady-State Free Precession (SSFP) MRI of the heart, cropped to focus on the heart, was warped with the φ_α's and altered with slight Gaussian noise (Figure 1).

Test Image	Whirl (α=50°)	Warped Image

Figure 1 Example of synthetic 3D images and whirl deformation fields represented by a warped grid

The eight warped images T were registered to the test image R using LogDemons and iLogDemons (σ_x=1mm, σ²=1mm and σ_f²=1mm). Relative mean square errors in grey level intensities (RMSE=||R-Toφ||_L2/||R-T||_L2), Jacobian determinant and distance to the true field φ_α (DTF = ||φ - φ_α||_L2) are reported in Figure 2. Deformation fields estimated by iLogDemons were almost incompressible. Jacobian determinants were always equal to 1±0.02 independently of the strength of the deformation to recover. Image matching accuracy was not affected by the incompressibility constraint (0.6% decrease). But most importantly, iLogDemons significantly improved the accuracy of the deformations. Mean and SD of DTF were systematically lower (average improvements of 29% and 36% respectively). The larger the deformation, the more significant the improvement. This points out the importance of the choice of the deformation model. Regions with homogeneous grey levels provided little information to accurately estimate the underlying deformation, as illustrated in Figure 3. The incompressibility constraint coped with this limitation by ensuring that the estimated deformation was of the same type as the true field. Similar conclusions were drawn on experiments with other parameters (σ_x=2, σ² = 2, σ_f² = {0.5,2}).

Figure 2 Results on images warped by 3D synthetic volume-preserving whirls with increasing whirl angle.

Figure 3 Streamlines of true and estimated whirl deformations (whirl angle α = 60°). iLogDemons provided incompressible deformations and outperformed LogDemons in terms of deformation field accuracy (yellow arrow)

References

1. Saddi, K.A., Chef d'hotel, C., Cheriet, F.: Large deformation registration of contrast-enhanced images with volume-preserving constraint. In: SPIE Medical Imaging. vol. 6512. (2007)
2. T. Mansi, X. Pennec, M. Sermesant, H. Delingette, and N. Ayache. LogDemons Revisited: Consistent Regularisation and Incompressibility Constraint for Soft Tissue Tracking in Medical Images. In Medical Image Computing and Computer Assisted Intervention (MICCAI), Lecture Notes in Computer Science. Springer, 2010. In press.

iLogDemons, Experiments on a Real Case

iLogDemons were then used to estimate the 3D left-ventricular myocardium strain from standard anatomical cine MRI of the heart. Such images have good in-plane and temporal resolutions but large slice thickness, which hampers the accurate estimation of cardiac through-plane motion (Figure 1). As the volume of the heart is almost constant during the cardiac cycle, incompressible registration is believed to improve the estimation of the cardiac deformation.

Figure 1 3D view of a standard anatomical cine MRI. Arrows represent the standard directions.
(Data courtesy of King's College London, St Thomas Hospital, Division of Imaging Sciences)

Anatomical short axis cine SSFP MRI (cMRI) of a patient with heart failure were acquired with multiple breath-holds (Achieva, Philips Medical System, 30 frames, 1.5mm² isotropic in-plane resolution, 10mm slice thickness). 3D tagged MRI (tMRI) were also acquired during the same exam (CSPAMM, 23 frames, 0.9mm isotropic resolution, tag size 3mm). No manual tracking of the tag grids was available since this task is extremely difficult due to the 3D motion. All the images fully covered both ventricles and no slice misalignments were detected. The cMRI were linearly resampled to get isotropic voxel sizes. The tMRI were spatially and temporally aligned to the cMRI using DICOM information. Because the transformations provided by demons algorithm are resampling fields, myocardium deformation was estimated by recursively registering all the frames of the cardiac sequence to the end-diastole (ED) time frame, as in (Mansi et al., 2009) (Figure 1). Registration parameters were σ_x=1mm, σ²=2mm, σ_f²=0.5mm (the smoothing was increased to accommodate the lower image quality). A 2-level multi-resolution scheme was used and registration was stopped as soon as RMSE stopped decreasing. Incompressibility constraint was applied only within the myocardium as volume of surrounding structures like blood pools vary.

Figure 2 Recursive tracking algorithm. Knowing the velocity v_In-1→I0 (in green): 1) Estimate v_In→In-1 (in blue). 2) Concatenate $v_In→In-1 and v_In-1→I0 (in grey). 3) Estimate v_In→I0 using 2) as initialisation (in red).

First, we estimated the myocardium motion by tracking the heart in the 3D tMRI using iLogDemons. For visual assessment, the deformations were applied to virtual planes manually positioned at ED (Figure 3). Realistic deformations consistent with the tag grids were obtained, which was further confirmed by the temporal variation of the radial, circumferiential and longitudinal myocardium strains (Figure 4, green curve). These results were similar to those obtained with logDemons, as the tag grids provided enough texture information in the myocardium to guide the registration. Hence, as no ground truth was available, we considered the iLogDemons estimation as reference. We then estimated the 3D motion of the heart from the cMRI and compared the results with the reference tMRI motion (Figure 4, blue and red curves). Visually, the warped virtual planes showed that incompressibility constraint did help to recover the longitudinal motion despite the large slice thickness of the cMRI. Estimated longitudinal and circumferential strains confirmed this finding. iLogDemons was closer to the reference than the logDemons (59% of improvement for radial strain, 84% for circumferential strain and 42% for longitudinal strain). Radial strain amplitude was more reasonable and the variations of the circumferential and longitudinal strains presented realistic patterns (Moore et al., 2000), where logDemons estimated an abnormal lengthening at the beginning of the cardiac contraction.

Figure 3 Close-up of the tMRI at end-systole with warped tag planes overlaid.
iLogDemons better estimates longitudinal and circumferential motion
(Data courtesy of King's College London, St Thomas Hospital, Division of Imaging Sciences)

Figure 4 Myocardium strains computed from short-axis cine MRI and tMRI. Mean and standard deviation computed over the entire left ventricle. iLogDemons better estimates longitudinal and circumferential strains.

Discussion and Conclusion

We have adapted logDemons algorithm to provide incompressible deformations. This has been possible by showing that demons Gaussian smoothing minimises an infinite order Tikhonov regulariser. This framework opens the way to new regularisers, such as elastic regularisation. As a result, incompressibility could be ensured by constraining the velocities to be divergence-free. The proposed incompressibility constraint does not introduce any new parameter. Those listed in this paper are present in any recent demons algorithm (Vercauteren et al., 2008). One could constrain the correspondence velocity to find the optimal incompressible update deformation. Yet, non-reported experiments showed that this does not significantly improve the results compared to iLogDemons: The updates are usually small and thus near-incompressible. The next step is to modify the demons energy to automatically handle incompressibility in subdomains of the image. From a clinical point of view, we are currently validating this method for the automatic estimation of 3D myocardium strain from standard cardiac images.

References

• T. Mansi, J.-M. Peyrat, M. Sermesant, H. Delingette, J. Blanc, Y. Boudjemline, and N. Ayache. Physically-Constrained Diffeomorphic Demons for the Estimation of 3D Myocardium Strain from Cine-MRI. In Proceedings of Functional Imaging and Modeling of the Heart 2009 (FIMH'09), volume 5528 of Lecture Notes in Computer Science, pages 201-210, 3-5 June 2009
• T. Mansi, X. Pennec, M. Sermesant, H. Delingette, and N. Ayache. LogDemons Revisited: Consistent Regularisation and Incompressibility Constraint for Soft Tissue Tracking in Medical Images. In Medical Image Computing and Computer Assisted Intervention (MICCAI), Lecture Notes in Computer Science. Springer, 2010. In press.
• C.C. Moore, C.H. Lugo-Olivieri, E.R. McVeigh, and E.A. Zerhouni. Three-dimensional systolic strain patterns in the normal human left ventricle: Characterization with tagged MR imaging. Radiology, 214(2):453-466, 2000.
• Tom Vercauteren, Xavier Pennec, Aymeric Perchant, and Nicholas Ayache. Symmetric log- domain diffeomorphic registration: A demons-based approach. In Medical Image Computing and Computer Assisted Intervention (MICCAI), volume 5241 of Lecture Notes in Computer Science, pages 754-761, New York, USA, September 2008. Springer-Verlag.

Inflammatory diseases

Summary

1. Clinical Introduction
2. Task Objectives
3. Data Collection
4. Multi-modal Image Registration [UPDATED]
5. Volume Estimation [UPDATED]
6. Erosion Assessment [UPDATED]
7. References [UPDATED]

Clinical Introduction

Juvenile Idiopathic Arthritis (JIA) is a very heterogeneous condition combining different forms of chronic arthritis of unknown origin. The onset of JIA ranges from several months after the birth to the age of 16, while the peak age is approximately at 6 years. JIA primarily affects the synovia and the bones of the joints; moreover chronic synovial inflammations may cause progressive joint destruction and serious functional disabilities. The impact of JIA on our society is significant, since it affects approximately one in 1000 children, and represents the major cause of acquired disability in paediatric population.

The causes of JIA have not been yet identified insofar that, in most cases, JIA is still a diagnosis of exclusion. Even if arthritis is not directly inherited, it seems very likely that both genetic and environmental factors are important for the manifestation of the disease. Hence, many efforts have been produced in the last few years in order to understand better the causes of the disease and, consequently, to introduce more effective (drug) treatments. More specifically, one major issue of interest is the possibility of distinguishing/classifying among different levels of the disease, and thus the ability of developing more specific and appropriate treatments for each subtype of patients. However, no statistically significant results have been obtained so far, and current classification - which is based on clinical criteria - is still unsatisfactory: many of the proposed subgroups don't appear homogeneous enough. It is still not possible, for instance, to distinguish effectively, early in the disease course, patients who are most likely to develop joint damage (and who therefore require a more aggressive treatment at an early stage) from those who are having a mild disease course. Finally, in clinical trials, drug efficacy may be judged only on clinical parameters, since physicians still lack quantitative measures that can allow for an early identification of both the progression of joint damage and the drug efficacy on disease progression.

From the considerations above, one can easily understand how important it is to identify reliable indicators of JIA activity, which can provide us with quantitative information about the disease. Actually, such goal is the main motivation of the research work outlined in this report, in which we exploit information contained in radiological images and we try to integrate it with clinical data.

Several imaging techniques have been tested on their ability to visualise and quantify both synovitis and erosions. Our hypothesis is that automatic analysis of the acquired images can lead to the desired quantitative indicators. However, not all the imaging modalities can be conveniently used for this purpose. In fact, for example, X-rays can show bone erosion, cartilage loss (indirectly, through joint space narrowing), and joint misalignment, but can not visualise synovia, joint effusion, articular cartilage, bone marrow, or ligaments and tendons directly. Thus, X-ray scoring methods can allow for joint damage evaluation in JIA, but are insensitive to synovitis.

Magnetic resonance imaging (MRI) is able to image synovitis and bone oedema/inflammation as well as damage to cartilage and bone, and it can detect erosive changes with greater sensitivity than radiography, particularly in early stages of the disease. MRI is also capable of detecting tendon pathology and evaluating ligament integrity.

Musculoskeletal ultrasound (MSUS) has been shown to be superior to clinical examination in the diagnosis and localisation of joint effusion, bursal fluid collection and synovitis. A number of studies have described improved sensitivity for detection of bone erosions in joints with the use of ultrasound as compared with conventional plain radiography. It is the imaging modality of choice for the diagnosis of tendon pathology. Ligament, muscle, peripheral nerve and cartilage pathology can also be readily investigated by MSUS. It allows a quick, safe and inexpensive access to otherwise undetectable anatomical information on the early targets of most rheumatic diseases.

MSUS and MRI are, therefore, potential powerful imaging tools to assess joint inflammation and the progression of joint damage. However, to date, little information exists on the use of these image techniques in JIA and standardised, validated, and feasible assessment systems are lacking.

The value of MRI, as an advanced method to evaluate disease activity and disease damage in adults with rheumatoid arthritis, is under active investigation by a research consortium called Outcome Measures in Rheumatology Clinical Trials (OMERACT). However, the results that will be drawn from OMERACT studies are not directly applicable to children, because adult rheumatoid arthritis is different from JIA and because the growing skeleton of children needs a different approach. Indeed, in children ossification is incomplete and joint space width varies with age.

Task Objectives

Data analysis in the scope of JIA pathology has two main objectives: 1) the identification of meaningful features that can be extracted from medical data, and consequently 2) the modelling of joint damage for an early diagnosis and assessment of the disease. As already pointed out above, our two main sources of information are clinical data and images. More specifically, the hypotheses we want to investigate is that a first set of feature can be extracted directly from radiological images of the joints, since we believe these are the most "characterising" feature (i.e. they discriminate better among patients). In a second stage of the classification process, image-based features can be integrated with those extracted from clinical data, in order to obtain a statistically grounded model of the disease.

Even if both image modalities, MR and MSUS, have shown to be very sensitive to early changes in the joints, our focus is only on the first modality, which appears better suited for quantitative measurements and for standardisation of the protocol. Henceforth, we assume MRI are the input images of our algorithms, and all the following discussions and assumptions apply specifically to this imaging modality.

From the wide literature on rheumatoid arthritis, we extracted the following set of features that have been reported to be effective in quantitative assessment of arthritis activity and damage.

• inflamed synovial membrane;
• bone erosion.

Therefore, from an algorithmic point of view, our main objective is to design, implement and validate algorithms to 1) (semi-)automatically extract and 2) quantitatively measure the above features from MR images.

The measurement of the synovial membrane volume is of particular interest. Indeed, such volume has been shown to be an early predictor for joint damage in adult patients affected by rheumatoid-arthritis (Østergaard et al., 1999) and osteoarthritis (Farrant et al., 2007), and it is expected to be valid also for JIA. The detection of bone erosions is important as well, since in MRIs, as opposed to conventional radiography, it can be successfully performed also in the early stages of the disease. Therefore both features, synovial volume and bone erosions, are important in modelling the disease course.

From a technical point of view, measurement of the synovial membrane volume completely depends on its segmentation. Working with 3D MRIs with no gap or overlap between the slices, as it is our case, the segmentation of the synovia from the other tissues allows the estimation of its volume by counting the number of detected voxels. The actual volume is given by the number of voxels times the volume of a single voxel. The quality of the estimate will depend on two factors: the accuracy of the segmentation, and the entity of the partial volume effect.

In the case of bone erosions we are interested in their detection, rather than in a quantitative estimate of their volume. However, the problem is in many respect similar to the segmentation of the synovial membrane, since it involves the classification of voxels or patches as belonging to an erosion or not. Differences will be in the features used to describe the voxels/patches, and possibly in the exploitation of methods to detect the bones contours (where the erosions are necessarily located).

Once one is equipped with the results of the two segmentation algorithms, it is possible to combine them in order to gain the confidence in the corresponding measurements. As it will be discussed in detail in the section about the data collection, for each enrolled patient there are MRI data at different modalities and different times (baseline and follow-up). Therefore, in order to combine the results, we also need to register the input images of the same patient.

To summarise, the achievement of our (clinical) goals implies the solution of the following (algorithmic) Computer Vision problems:

• synovial membrane segmentation in contrast-enhanced MRIs;
• detection of bone erosion in T1 pre-contrast MRIs;
• registration of intra-patient images at different times and different modalities.

In the following subsections we provide brief overviews of related works for each of the technical objectives.

Estimation of the Synovial Membrane Volume

The segmentation of the synovial membrane has often been performed manually by human operators (Østergaard et al., 1999; Farrant et al., 2007). However such approach is extremely time-consuming, and furthermore the reliability of the results can suffer from intra- and inter-observer variability. Even if it is widely acknowledged that automatic segmentation can be advantageous in terms of both time and consistency of the results, relatively few works propose automatic methods for the segmentation of the synovial membrane, especially for the sites of our interest, hips and wrist.

The automatic segmentation problem can be approached from two different directions: either by looking for the contours/surfaces of the regions in 2D/3D image (i.e. by contour-based segmentation), or by labelling directly the pixels/voxels (i.e. by pixel-based segmentation). By surveying the tiny state-of-the-art works on the segmentation of the synovial membrane, the pixel-based approach results the leading one. Surprisingly enough, so far the focus has been only on the simplest pixel-based segmentation method: the thresholding in its various forms. In global thresholding the intensity range of the images is split into intervals and each pixel is labelled according to the interval its intensity belongs to. A comparison of segmentation of the synovial membrane by thresholding with manual annotation is reported in (Østergaard, 1997). This technique is obviously too basic, and necessarily is highly sensitive to noise and artifacts in the images, inhomogeneities in the magnetic fields. Moreover if one looks at the value of a single pixel alone, the spatial coherence between neighbouring pixels cannot be taken into account. Another widespread segmentation method in medical image analysis is fuzzy c-means clustering, and it has been recently applied to the segmentation of the synovia in the wrist (Tripoliti et al., 2007).

For completeness, we now report some recent works on the segmentation of the cartilage, because this problem shares some affinities with our. (Folkesson et al., 2005a) and (Folkesson et al., 2005b) segment the cartilage in MR knee scans, with the aim of quantifying the progression of osteoarthritis. The proposed methods are both based on the combination of Approximate Nearest Neighbour (ANN) classifiers; an ANN is similar to a plain k-Nearest-Neighbour classifier, but it looks for an approximate solution in order to reduce the computational cost. The method presented in (Grau et al., 2004) has also been applied to knee cartilage. It extends the popular watershed segmentation algorithm by introducing prior knowledge obtained by estimating the posterior probabilities of the classes under study assuming a Potts' model, a generalisation of Ising model for two classes. (Warfield et al., 2000) presented a semi-automatic segmentation method, applicable to a wide range of tissues (among them also cartilage), iterating between a classification step and a template registration step. At the classification step each voxel is assigned to the class most represented by its first k neighbours in the aligned template (which has been previously segmented). The intermediate classification result is then used to refine the non-rigid registration. (Tang et al., 2004) applied a modified version of the gradient vector flow (GVF) snakes to the segmentation of ankle cartilage. In (Solloway et al., 1997) the active shape model is applied, unmodified, to the segmentation of femoral cartilage.

Detection of Bone Erosions

In order to detect bone erosions, two different approaches have been proposed so far. First, since it seems plausible to consider the erosions as “lesions” (or abnormalities) connected to bone boundaries, a two step procedure can be adopted in which the entire bone volume is pre-segmented, and then the abnormalities are searched along the boundaries of the segmented volume. The alternative approach aims to detect directly the regions in the images where bone lesions are most likely to occur. The first approach has drawn more attention than the second one. Indeed, most of the previous works focused on finding the contours of the bones in plain images (i.e. each MR slice separately). However, they share the disadvantage that the performance of the detection is dependent on the accuracy of the delineation of the contours. In order to increase such accuracy - and, thus, to reduce the risk of misclassification - the algorithms usually rely on human initialisation and intervention.

On the other side, a direct detection of erosions necessarily requires both a formal definition of what an erosion is - i.e. how it looks like in an MR image - and some assumptions about where it can happen; otherwise the number of false detections can make the performance of the system fall down. This is because the local appearance of bone erosions is indistinguishable from many other artifacts that do lie within the wrist area. As for the other approach, the human contribution may improve significantly the performances of the system.

Therefore, if one aims at comparing and evaluating different methods for bone erosion detection, one important aspect that must be taken into account is the degree of autonomy, which refers to the need for human intervention either before or during the running of the algorithms. Intuitively, the physicians might be asked to assist the algorithm only rarely and if no other options are available. In fact, human intervention is shown to be critical when physicians are required to manually segment image regions, due to the great variability of their responses. Indeed, it is still matter of debate how to combine different manual segmentations performed by either the same physician several times or two or more physicians. Bias can be introduced if this step is not performed carefully. Given all these caveats, we strongly believe that human intervention should be avoided almost completely. Therefore, later in this report we will propose a method to overcome this problem, by making the segmentation process more complicated but, at the same time, more robust.

The following is a more detailed discussion of some of the most popular algorithms that have been used so far to detect bone erosions.

Let's consider first the methods based on active contours [Kass et al. 1988]. In a nutshell, active contours look for a 2D curve minimising an energy made up of two terms: an internal energy which accounts for the elasticity and smoothness of the curve, and an external energy which accounts for the image features (for instance the edges extracted by Gaussian gradient filtering). Note that although the solution is highly dependent on the model initialisation, the method might be made more robust by adopting a coarse-to-fine strategy.

This general approach has been applied, with minor modifications, to the segmentation of bones in both radiographies and MRIs. The former is the case of the work described in [De Luis-Garcia et al. 2003], where the authors require an accurate segmentation of the finger bones in order to estimate the skeletal age of the patient. Their method modify the general approach in two ways: first, they exploit domain knowledge to initially place the snakes on the phalanges and metacarpal bones; second, the external energy is the sum of the image gradient and of an energy which favours the expansion of the contour (introduced in [Cohen 1999], to make the method more robust to initialisation). Active contours have been applied to bone MRIs in [Snel et al. 1998], who looked for the planar contours of the carpal bones in the MRI slices. Inspired by the work described in [Lobregt et al. 1995], they use only the normal component (with respect to the contour) of the external energy, which avoids the problem of clustering the pixels of the curve (henceforth called snaxel), which affecting the original method. Moreover, in order to overcome the problems due to the varying image contrast along the contour, they also normalise the external energy with a local contrast level estimated in the neighbourhood of each snaxel. In both works, the snake can be refined or snaxels fused together depending on the distances between each of them.

Active contour methods do not use any prior knowledge on the shapes of the contours. Hence a more sophisticated class of methods has been introduced, called Active Shape Models (ASM), which try to exploit any prior information available for the specific problem, by analysing the mode of variations of a set of representative examples (the so-called training set). In the Medical Imaging context, ASMs have been already applied to the segmentation of vertebrae in Dual Energy X-ray Absorptiometry (DXA) images [Smyth et al. 1997] to overcome the problem posed by the low signal-to-noise ratio and poor spatial resolution of DXA with respect to standard radiographies. The authors claim an accuracy comparable with that of human operators. In [Behiels et al. 2002], an assessment is presented, on 400 images of different bone structures, of the performance of three different methods: a standard ASM method, a probabilistic method based on Maximum-a-Posteriori criterion [Wang et al. 1998] and an extension of ASM. This third method is essentially based on the introduction of a regularisation term in the fitting problem, which is then solved by minimal cost path (MCP) search.

Registration of multi-modality MR images

As we made it clear above, the results obtained separately by the two segmentation algorithms (synovial membrane and bone erosions) would improve if one was able to compare (spatially) the locations of the segmented regions. In fact, inflammation of synovia seems to be caused by specific neurotransmitters released from inside the bones in the regions where an erosion occurred. However, since our aim is to combine statistically meaningful (quantitative) features extracted from heterogeneous data sources (clinical and images), it is possible that the advantages of working on registered images, rather than on the measures extracted from unregistered images, can be quite limited. In fact, other clinical factors could carry similar information insofar that information related to the registration would result in more redundancy of the dataset. Therefore, so far we considered the registration of multi-modality MR images more as an exploratory tool than as an essential component of our statistical decision system.

The actual work we have done so far focuses on surveying the scientific literature on multi-modality registration of medical images. A complete report on the state-of-the-art works is still in preparation and will be available in next documents of the project.

Data Collection

Planned Data Acquisitions

The Health-e-Child study is designed as cross-sectional and prospective, with an expected patient sample of 300 JIA patients (100 in each of the 3 participating centres). Patients with JIA (according to ILAR classification) and with active arthritis in the wrist and/or hip are enrolled in the study. For enrolled patients, MRI and US of the wrist and/or hip are performed at study entry (baseline), and then at one and, when possible, two years of follow-up. MRI of the wrist is performed only in cooperating patients that do not require general anaesthesia.

The images are acquired with the protocol developed in the first year of the project, consisting of:

• volumetric, fast field echo T1
• volumetric, fast field echo T1, with fat suppression, 1.5 min after injection of Gd-DTPA
• turbo spin echo T2, with fat suppression
• volumetric, fast field echo T1, with fat suppression, 10 min after injection of Gd-DTPA

Fig.1 - Examples of the different types of acquisitions from the MR protocol. The images in the same column are slices at the same position along the Y axis.

	y=-70	y=-65	y=-60
T1 pre-contrast
T1 after 1.5 min
T2
T1 after 10 min

As explained in the section about the task objectives, the results of the image analysis will be integrated with other data in order to provide a model of the disease course. In this early stage of the research work, we are mostly interested in integration with clinical parameters (including standardised measures of disease activity and disease damage), routine laboratory data (including ESR, CRP, WBC and differential count, RBC, Hb, PTL, ANA, RF), and semi-quantitative scores for damage evaluation in XR, US and MR.

Data Annotations

Our approach to image analysis is based on machine learning; more specifically we design our algorithms in the setting of learning by example. Therefore, we need annotated examples, highlighting the structures of interest, in order to perform training, validation and testing of the methods we develop. Namely, we need a manual segmentation of the synovial membrane for the volume estimation and of the bone erosion for erosion detection. As we will see in short, the former annotations must be particularly accurate in order to obtain an accurate estimation of the volume. For the annotation task we developed a software tool described in the deliverable for WP12.

Fig.2 - Examples of manual annotations from two operators. The images in the same column are slices at the same position along the Y axis.

	y=-70	y=-65	y=-60
Operator 1
Operator 2

The analysis of the manual annotations performed on the same exams by two clinicians provides useful indications on the complexity of the problem. First, we can compare the annotations on a per-slice basis. If we compare the synovial membrane volume resulting from the annotations, we can observe that: 1) they are highly correlated, with r=0.92 and p<0.001, but 2) there is an average absolute difference of 30.5%. The latter measurement is mitigated when we look at the median absolute difference, which is 14.4%, but the fact remains that the volumes estimated from different observers shows a variability. This variability does not disappear considering whole exams: some of the differences tend to balance, and the absolute average difference drops to 15.2%, but the median absolute difference decreases only to 12.1%. Also interesting is the comparison between annotations from a classification point of view. If we look at the voxels which have been classified by both observers as belonging to the synovial membrane, these make up 63% (on average) of the voxels marked as synovia by either of them. That is, around one third of the voxels classified by one observer as synovia, are classified by another one as NOT synovia.

This results, although partially due to operator errors, are also the product of an inherently difficult problem. The structures to be segmented are extremely thin structures (see Fig.2), and often can be discriminated only on the basis of anatomical knowledge rather than signal intensity. As a result of the thin shape of the structures, small pixel-wise errors on their boundaries might produce significant errors on the volume. This can be verified by subjecting the manual annotations to binary dilation and erosion, and comparing the resulting volumes with the original ones. In both cases there is a absolute change of over 30%.

The analysis of the manual annotations, briefly summarised above, points out the complexity of the problem of accurately estimating the synovial membrane volume. It also sets out what might be a reasonable goal for the performance of an automatic method, which of course cannot beat the examples it is trained with.

Multi-modal Image Registration

Both clinical objectives of synovial volume estimation and bone erosion detection would benefit from the preliminary registration of the multi-modal MR images included in each study. In particular, we face three types of registration problems:

• intra-patient, intra-study registration
• intra-patient, inter-study registration (between baseline and followup studies)
• inter-patient registration

The first problem is the simple registration of all modalities in a single MR study. This problem can be approached as a rigid registration problem, and its solution is needed in order to be capable of fully exploit all information provided by the different modalities. Indeed, this type of registration allows to locate the image voxels corresponding to the same physical point across modalities. The second problem occurs when considering pairs of baseline and follow-up studies. In this case it would be helpful to be able to find the voxels corresponding to the same anatomical point, in order to perform visual or automatic comparison of image features indicating, eg. bone erosions. This problem is considerably more difficult than the previous one, since the natural growth together with the pathological processes produce deformations between the two studies which cannot be modeled as rigid. Finally, addressing the third problem would provide a standard frame of reference for all studies, which could potentially allow the exploitation of information about the voxels position both in synovia segmentation and in erosion detection. Moreover, a non-rigid inter-patient registration could allow for the use of an anatomical prior defined on a reference atlas, thereby improving the robustness of the segmentation.

All problems have been addressed. We employ a registration algorithm based on the maximization of the mutual information, as published in Mattes et al, 2003. The implementation is based on the open source library ITK (see http://www.itk.org). Note that although technically in the second and third problems the modalities are the same, in practice we found that maximizing the mutual information provides an efficient way to be robust with respect to the changes in image intensities and anatomy naturally occurring in this case. We show below some representative examples of the results obtained on our data for the intra-study registration and the inter-patient registration (only the rigid part). It can be readily seen that the three modalities are perfectly aligned. The inter-patient alignment also worked well, keeping in mind that pure rigid registration cannot compensate for the differences inherent to the different anatomies between patients.

Fig.3 - Examples of registered slices, at the same position along the z axis.

	modality #1	modality #2	modality #3
Exam #1
Exam #2

In the examples above only the a rigid transformation has been used. To address the problem of inter-study registration, in particular having in mind the application to the assessment of erosion progression, a non-rigid component is also required. The same registration method cited above, through the maximimization of mutual information, proved extremely effective. We show below two representative examples of registration between baseline and follow-up studies of the same patient, as well as between studies of different patients.

Fig.4 - Examples of non-rigid registrations. Slices at same position along z.

	study #1 (fixed)	study #2 (registered)	study #2 (original)
inter-study
inter-patient

We remark that, depending on the application, we use the deformation field estimated during the registration process, rather than the image warped with such field. This will be seen more in detail in the section about erosions.

Volume Estimation

The accurate estimation of the volume of the inflamed synovial membrane is the first objective of our work. From an algorithmic point of view, such objective can be achieved by first segmenting inflamed synovia from all others tissues, and then evaluating the volume of the segmented regions. The conversion between the volume in voxels of the segmented regions and their actual volume in cubic mm can be conveniently done by means of MR volumetric acquisitions. To this aim, the clinical partners designed a protocol for MRI acquisitions, including series of T1-weighted, volumetric MR images, of the hips or wrists of the enrolled patients. Two of the series are acquired 1.5 and 10 minutes after the injection of a Gd-DTPA bolus, which enhances the contrast of highly vascularized tissues, as it should occur in case of an active synovial pannus.

As we have seen in the previous section, the manual segmentation of the enhanced parts is, at least in the wrist, a tedious and time-consuming, and the results have an high inter-operator variability. These are serious hurdles to the use of the synovial volume as a predictive marker for the future erosive damage. Therefore, any automatic or semi-automatic measurement method would be welcome, provided it could reduce both the variability and the time required for the measure. However, one should not overlook the difficulties in designing such a method. These are due to the thin shapes under study (high sensitivity of the result to errors in the segmentation), and to the necessity of incorporate some prior anatomical knowledge (for discriminating between synovia and other enhanced structures).

We approach the problem as one of classification. Previously, we approached the problem in two stages. On a basic level, we wanted to classify the voxels as belonging to enhanced structures or not. Note that at this stage we were not yet interested in the discrimination between actual synovial membrane and other enhanced structures. At a second stage, however, once the enhanced voxels had been classified/segmented, we wanted to classify each group of connected voxels as belonging to the synovial membrane or not. This approach has been now superseded. Both levels of classification (enhanced vs. not enhanced, inflamed synovia vs. other enhanced structures) are addressed by a single classifier.

The classifier is obtained by training a discriminative function, represented by a linear combination of features selected from a large data-dependent dictionary. The training step is performed by solving a regularised least-square problem, while the feature selection is currently performed via Orthogonal Matching Pursuit (OMP, see Davis et al, 1997).

The performance of the classifiers and of the volume estimation is tested in a cross-validation fashion. We iteratively put aside one or more annotated exams and train/validate the classifiers on the remaining data; the obtained classifiers are then tested on the untainted exams.

Further details can be found in [Basso et al, 2009a], [Basso et al, 2009b] and [Basso et al, 2010].

Fig.5 - Examples of segmentation. The images in the same column are slices at the same position along the Z axis.

Patient #1
Patient #2 (baseline)
Patient #2 (follow-up)

Previous studies estimated for the volume estimation a precision of 0.025 cc (computed per slice, over the cross-validation loops). For comparison, a human observer compared with a second one would achieve a precision of 0.017 cc. In order to evaluate the usability of the system in a clinical setting, preliminary to validation, we computed its agreement with respect to the volume computed by a clinician. The system achieves an ICC rate of 0.85 (95% CI: 0.83-0.87), compared to the agreement between two human observers of 0.94 (95% CI: 0.92-0.95). This two results are excellent, especially since they have been obtained from a completely automatic system, without human intervention.

Further cross-sectional and longitudinal studies confirmed the viability of the synovial volume normalized by the body surface (NSV), as an imaging biomarker for assessing the disease activity and in studying its progression. A significant difference in the values of the NSV, from 3.48 to 1.79 (p<0.0005), has been found between baseline and follow-up, with a standardized response mean (SRM) of 1.04. A strong and significant correlation is present between the radiographic damage progression in patients undergoing DMARD therapy and the volume at baseline (r=--0.82, p<0.02). The cross-sectional analysis of 58 studies showed significantly different values of the NSV in patients grouped according to: low/high score pain (2.26/3.92, p<0.0002), low/high score swelling (2.45/4.28, p<0.0005), low/high score limitation of motion (LOM, 2.42/3.61, p<0.01), as well as polyarticular/non-polyarticular JIA subtype (2.33/3.30, p<0.02). Moreover, the NSV correlates (r=0.39, p<0.005) with a semi-quantitative score of synovitis based on the adult OMERACT score, which assesses both the thickness of the synovial membrane and its hyper-intensity. The NSV presents also moderate/weak significant correlations with all the clinical indicators of disease activity that the OMERACTsynovitis score correlates with: the physician score on a visual analog scale (r=0.42, p<0.002), the number of joints with limitation of motion (r=0.37, p<0.005), the number of swollen joints (r=0.36, p<0.01), the composite score JADAS71 (r=0.35, p<0.01), the white blood cells count (r=0.35, p<0.01).

Erosion Assessment

As anticipated above, bone erosions play a fundamental role in the analysis and quantification of JIA progression. Specifically, an accurate quantification of temporal changes of erosions between baseline and followup studies may be used to discriminate the actual outcome of specific drug treatments or therapies.

From a radiological viewpoint, an erosion may be defined as a bony depression visible in at least 2 planes on 3DT1 sequences. However, this basic definition may be equivocal, and usually radiologists need additional findings suggestive of inflammation in order to assess correctly the bone damage. Examples of such further findings are the presence of surrounding bone edema, effusion and enhancement in the corresponding areas on T2 sequences or the presence of irregularities on X-ray examinations). Most of such findings require an high-level analysis and the combined visual inspection of different examinations, which are very difficult to replicate using a computerized system. Nonetheless, despite being considered more achievable, also the manual analysis has shown to be affected by a large variability when performed by different physicians.

Fig.6 - An T1 pre-contrast image of a wrist (left). The red square represents the region of interest for the bone erosion modules, and is magnified on the right.

In collaboration with the physicians involved in the project, we considered the problem of defining a suitable and reliable scoring system for MR studies. It is worth emphasizing here that the main goal of this part of our work was to assist the physicians to improve the current diagnostic analysis of JIA patients. Therefore, the starting point has been the evaluation of the existing scoring system based on a semi-quantitative index ranging from 0 to 4, that shows a good reliability. However, the limited number of possible values makes the index poorly sensitive to changes, especially in children.

Fig.7 - Progression of the erosion after one year. The first two images (on the left) represent the baseline examination. The third image (on the right) show the followup.

The first result obtained within this part of the project was to propose a new score from 0 (normal) to 10 (the whole bone eroded). Of course, such score requires a more accurate analysis since the increased number of classes may lead to a poor reliability. The automation of the analysis in this context is therefore essential. In the last months of the project, within WP11, our work focused on designing, implementing and testing algorithms for the assessment of the new score and its monitoring over time. Specifically, we first developed a software suite quiped with an intuitive graphical interface enabling the physicians to exploit the results of the registration for a semi-quantitative evaluation of erosion progression. Thereafter, we worked on automating such evaluation by leveraging on the annotations collected by the physicians. The first part of our work comprised two distinct modules: one was related to the registration of baseline-followup pairs of images acquired on the same patient, while the other consisted in developing a user interface starting from the definition of suitable use cases taking into account the requirements of the physicians.

A considerable amount of our work has been done to adapt and exploit the results of task T11.3.3. Details about the registration part are discussed above. The accuracy of the overall system (both during the manual and the automatic analysis) strongly depends on the results of the registration. Indeed, by using the registered images (instead of the original baseline-followup pairs) we observed an increase in the reliability of the new scoring method.

Once we registered the baseline-followup images of all the patients, the second important step was to develop a system by which the radiologists may exploit the registration results. It is well known that the non-rigid registration modifies the anatomy of the tissues and, specifically in the case of bones, it may prevent a correct evaluation of morphological abnormalities. Having this problem in mind, we decided to introduce a navigation module in our system to allow the user to move from a baseline to the corresponding followup 3D study by simply pointing on a specific voxel. The following figure represents a screen-shot of the designed GUI.

Fig.8 - A screen-shot of the GUI used for the annotation of MR images of wrists in order to collect positive and negative examples for the classification module.

When a user clicks on a voxel (e.g. one corresponding to a suspect erosion) on the baseline 3D image (the leftmost image panel in the GUI) the system automatically updates the panel on the right to display the corresponding slide in the followup image. Furthermore there will show up a small panel displaying a zoomed neighborhood of the selected pixel in both baseline and followup images. Such visualization scheme is useful because it exploit the results of the registration without modifying the structures of interest in the image.

After an initial stage of testing by the clinicians, the tool has been improved by allowing the users to annotate the exam directly through the graphical interface. This led to a reduction of the time needed for fully annotating a single pair of exams. Furthermore, this gave us the opportunity of collecting data for the subsequent tasks, i.e. the further automation of the analysis.

In order to accomplish the automatic analysis of the bone erosions using machine learning techniques, one needs to collect a large training set of both positive and negative annotated examples highlighting the structures of interest. In the context of JIA, annotated examples may be voxels in the 3D MR volumes of T1 pre-contrast wrists. Positive voxels are those in the neighborhood of an erosion while all the others in the volume are negatives. For collecting the training examples we implemented some slight modifications to the software in order to save the annotations during the ordinary analysis. Despite the ease of saving the annotations, the size of the training set is not large because the number of patient enrolled in the clinical study was limited.

Our approach to the automatic analysis of erosion is based on a two stage schema, in which we first identify some coarse regions of interest (connected the contours of bones) in the MR sequence and then detects those group of voxels that are likely to represent an erosion. The first segmentation is required to reduce the number of false positives and to reduce the computation time of the whole algorithm by discarding part of the (very unlikely) candidate erosions.

In order to accomplish the classification step, we decided to adopt the same learning method designed for the synovia, and the feature extraction algorithms have been used to discover the optimal set of cues and features for this specific task. The set of higher-level features are based on the following cues:

• Intensities of the voxel and the set of its first neighbors.
• Position of the voxel in millimeters, relative to the center of the acquisition volume.
• Multiscale 2-jets, computed by convolution with the appropriate derivative of a Gaussian kernel.

Fig.9 - An example of a slice with annotated erosions.

References

• C. Basso et al, 2009a. Segmentation of Inflamed Synovia in Multi-Modal MRI. In Proc. of IEEE ISBI 2009, June 28 - July 1 2009.
• C. Basso et al, 2009b. Automatic Annotation of 3D Multi-Modal MR Images on a Desktop Grid. MICCAI-Grid Workshop, held at MICCAI 2009, 20-24 September 2009, London, UK.
• C. Basso et al, 2010. Automatic estimation of inflamed synovial membrane volume in 3D MR images. Poster at ECR 2010 Scientific Exhibit, 4-8 March 2010, Vienna, Austria.
• C. M. Bishop, 2006. Pattern recognition and machine learning. Springer Series on Information Science and Statistics. Springer, New York, NY, USA.
• Davis et al., 1997. Adaptive greedy approximations. Constr Approx, 13(1):57-98.
• J. Farrant et al., 2007. Advanced imaging in rheumatoid arthritis. part 1: Synovitis. Skeletal Radiology, 36(4):269–279.
• J. Folkesson et al., 2005a. Combining Binary Classifiers for Automatic Cartilage Segmentation in Knee MRI. Proc. CVBIA 2005, LNCS 3765, pp. 230-239.
• J. Folkesson et al., 2005b. Automatic segmentation of the articular cartilage in knee MRI using a hierarchical multi-class classification scheme. Proc. MICCAI 2005, pp. 327-334.
• V. Grau et al., 2004. Improved watershed transform for medical image segmentation using prior information. IEEE Trans. Medical Imaging, 23(4):447-458.
• L. Lo Gerfo et al., 2007. Spectral algorithms for supervised learning. Submitted to Neural Computation. Preprint available from rosasco@disi.unige.it
• Mattes et al., 2003. PET-CT image registration in the chest using free-form deformations. IEEE Trans. Med. Imaging, 22(1):120-128.
• M. Østergaard, 1997. Different approaches to synovial membrane volume determination by magnetic resonance imaging: manual versus automated segmentation. British Journal of Rheumatology, 36(11):1166-1177.
• M. Østergaard et al., 1999. Magnetic resonance imaging-determined synovial membrane volume as a marker of disease activity and a predictor of progressive joint destruction in the wrists of patients with rheumatoid arthritis. Arthritis & Rheumatism, 42(5):918–929.
• S. Solloway et al., 1997. The use of active shape models for making thickness measurements of articular cartilage from MR images. Magnetic Resonance in Medicine, 37(6):943-952.
• J. Tang et al., 2004. Ankle cartilage surface segmentation using directional gradient vector flow snakes. Proc. of ICIP'04, pp. 2745-2748.
• J.P.Thirion, 1998. Image matching as a diffusion process: an analogy withMaxwell's demons. Medical Image Analysis, 2(3):243-60.
• E. E. Tripoliti et al., 2007. Automated segmentation and quantification of inflammatory tissue of the hand in rheumatoid arthritis patients using magnetic resonance imaging data. Artificial Intelligence in Medicine, 40:65-85.
• S. K. Warfield et al., 2000. Adaptive, template moderated, spatially varying statistical classification. Med. Image Anal., 4(1):43-55.

Brain tumors

Integrated disease Models for Paediatric Brain Tumours

I- Tumour Growth Modelling in the Brain

1. Literature Review of Tumour Growth Models
2. A Macroscopic Patient Specific Growth Model
3. An Eikonal Approximation for the Reaction-Diffusion Type Tumor Growth Models
4. Extrapolating Invisible Tumour Distribution
5. Parameter Identification Problem
6. Monitoring Slowly Evolving Tumors

II- Computational Anatomy of the Brain

1. Computational Anatomy: Aims and Methods
2. Morphometry of the Cortex Inferred from Sulcal Lines
3. A Statistical Framework for Anatomical Curves and Surfaces
4. A Generative Model of Brain Populations Variability
5. NEW Spatiotemporal Atlas Estimation in Longitudinal Datasets NEW

Tumour Growth Modelling in the Brain

1. Literature review and general overview

1.1 Introduction

Tumour growth modelling is the study of complex dynamics of cancer progression using mathematical descriptions. Internal dynamics of cancerous cells, their interactions with each other and with healthy tissue, nutrition and oxygen transport from the extracellular matrix (ECM) and from the vascular network, chemicals secreted by tumour cells, type of the underlying tissue and many other factors are modelled by mathematical abstractions. These abstractions rely on biological and clinical observations coming from different sources at different scales. Example observations include in-vitro experiments, in-vivo experiments on animal subjects, biopsy and autopsies results, medical images of patients like computed tomography (CT) scans or magnetic resonance images (MRIs). These experimental results and images are keys to develop models describing the tumour growth process accurately.

Mathematical growth models offer important tools for both clinical and research communities in oncology. By providing a common mathematical ground to combine experimental findings made in many diverse fields of cancer research they allow us to interpret such experimental results and understand the underlying mechanisms of tumour growth. Virtual experiments and simulations give the opportunity to observe effects of different treatments on cancerous cells, and could lead to improve these treatments or suggest new ones. Personalised models, adapted to patient specific cases, could be used in therapy planning by suggesting irradiation regions adapted to growth dynamics or optimal temporal distribution of chemotherapy.

1.2 Mathematical Growth Models

Tumour growth modelling has received considerable attention in the literature and there has been a vast amount of work during the last 15 years. A complete review of all the work done on this subject is outside the scope of this page. Here, we synthesise the main approaches towards this problem and review some of the latest work, in view of giving the reader an overview of the current trends.

Tumour growth is an extremely complicated process consisting of different interacting sub-processes at different scales. Proposed models concentrate on one or several of these sub-processes. Research conducted on tumour growth modelling can be coarsely classified into two large groups: microscopic models and macroscopic models. The main difference between these two classes is the scale of observations they are trying to explain and formulate. Microscopic models concentrate on observations in the microscopic scale, like in-vitro and in-vivo experiments. As a result they formulate the growth phenomena at this scale. Macroscopic models on the other hand, concentrate on observations at the macroscopic scale like the ones provided by medical images. They formulate the average behaviour of tumour cells and their interactions with underlying tissue structures, which are visible at this scale of observation (grey matter, white matter, bones...). These models try to describe the behaviour of the tumour as a whole, consisting of millions of cells.

Further classification within these groups can be made based on the stage of the tumour growth being analysed (avascular growth, angiogenesis and vascular growth) or the effect of the growth on the brain (diffusion/invasion of tumour cells and mass-effect). In this review we use the classification based on the stage for the microscopic models and the one based on the effect for the macroscopic models for better comprehension.

Classification of tumour growth models: First level classification is based on the scale of the model and scale of the observation model uses. Second level classification is based on the stage of the tumour growth for microscopic models and effect of the tumour for the macroscopic models.

1.2.1 Microscopic Models

Microscopic growth models aim to describe the tumour growth process at the cellular level based on the experimental observations at this level. They take into account physical and chemical interactions between cancer cells, the extra-cellular matrix and the healthy cells. Mechanical phenomenon like pressure, cohesion and adhesion forces are often included in these descriptions. As for chemical interactions, microscopic models formulate phenomena like diffusion of nutrition and oxygen, secretion of different factors by tumour cells and their effects on the surrounding. Mathematical systems obtained for these models are usually very detailed because they try to take into account all the observed factors. As a consequence, the number of parameters for such models are very high. From the technical point of view, formulations used in microscopic models enjoy a large variety of mathematical methods. Most commonly used methods are partial differential equation (PDE) systems, cellular automata and statistical models.

• Avascular Growth/Solid Tumour:

The avascular growth corresponds to the stage where the growth process is mostly governed by the proliferation of tumour cells. In this stage the tumour is considered to be a solid mass, which is growing by mitosis. Although not completely known, it is thought that there is no invasion of the healthy tissue. The interactions between tumour cells and the healthy tissue is also thought to be limited [Araujo and McElwain, 2004]. The tumour cannot grow indefinitely in the avascular stage because as the tumour mass grows, less and less nutrition is available for the cells deep inside the avascular mass. As a result, necrosis begins: tumour cells that are not getting enough nutrition die, and only cells on the outer perimeter of the tumour continue to proliferate. At one point necrosis and the proliferation balances each other and the avascular tumour reaches a limit size, which is assumed to be around 1-3 mm in diameter [Orme and Chaplain, 1996].

Most growth models proposed to formulate this stage try to describe mathematically three different experimental observations: the formation of the necrotic core, dormancy of the tumour at the limit size and the motion of cancerous cells within the solid tumour. In order to achieve this, they formulate different causes like:

• Mitosis/Population growth,
• Distribution and consumption of oxygen,
• Apoptosis/Programmed cell death,
• Surface tension,
• Cell cycles,
• Secretion and diffusion of growth inhibitory hormones (GIFs),
• Chemotaxis/Active migration,
• ...
Building blocks of these models are coupled ordinary and partial differential equations such as population growth and reaction-diffusion. Variation of parameters of these equations have also been studied to assess the importance of deviation from the mean behaviour.

• Tumour-Induced Angiogenesis:

This is the stage where tumour cells in the avascular mass modify the existing vascular structure to create new vessels that would feed them. Through this process the tumour can overcome its limit size, grow much faster and invade the surrounding tissue. Due to this crucial role of angiogenesis, the underlying mechanisms of this stage have captured attention and many models have been proposed.

Tumour-induced angiogenesis is a very complex process including lots of chemical and mechanical phenomenon, which have not been totally understood. [Mantzaris et al., 2004] reviewed some of the known biological processes taking place in angiogenesis. The basic observable consequence is that tumour cells affect nearby blood vessels to sprout new vessels towards itself, creating new vascular structure. As a result of angiogenesis, the tumour receives extra nutrients and oxygen and the growth speeds up. Mathematical formulations of this stage usually takes into account the motion of endothelial cells (ECs), tumour angiogenesis factor (TAF), fibronectin structures of extracellular matrix (ECM), vascular endothelial growth factor (VEGF) and angiogenesis inhibitor factors.

As in the case of solid tumour growth, many angiogenesis models are based on integro-differential equations. In most of them the diffusion of the above-mentioned chemicals and the density of blood vessels are formulated using reaction-diffusion equations. In addition to these, some discrete models are trying to predict the structure of the formed vasculature in 2D, rather than formulating it as a density. Such discrete models use different tools like Nash games and stochastic differential equations.

• Vascular Growth/ Invasive Tumour:

The third stage of the tumour growth, vascular growth, has been paid less attention than the two previous stages. The complexity of the tumour growth in this stage is higher because there are several processes going on simultaneously. In addition to cellular and chemical interactions going on as in first two stages, tumour cells start to invade the surrounding tissue via mechanisms that are not clearly known yet. At this stage, the tumour becomes diffusive and is not considered to be solid anymore. While the difference between cancerous and healthy regions are clear in the avascular stage, this difference vanishes during the vascular growth because tumour cells move towards healthy regions. Vessels inside the vascular tumour might have been formed by angiogenesis or the tumour might have initialised around a vessel, as in the case of tumour cords. Unlike avascular tumours, the source of nutrition is not limited to diffusion from the perimeter. Thus, the formation of necrotic regions is much more complex, if they exist at all. For the same reasons, vascular tumours are not compact masses of cancerous cells, they don't have a limit size and they can grow indefinitely.

First models examining the vascular invasion concentrated on the metastases process (invasion through blood stream). Using similar ideas, the local invasion is modelled by formulating the infiltration of the tumour cells. Different models have taken different approaches for this. Some groups modelled the infiltration as directed diffusion, particularly towards blood vessels, or as an active transport, resulting in a reaction-diffusion formalism. At the opposite of these continuum models, some cellular automata models were also proposed for the dynamics of cancerous cells, along with a reaction-diffusion formalism for the chemicals taking a role in the growth process. Besides the methodological differences with the previous models, new phenomena like the effect of the acidic metabolism of tumour cells or the effect of the vascular structure on the tumour growth are studied and included in vascular/invasion models.

1.2.2 Macroscopic Models

Observations at the macroscopic scale consists of medical images like computed tomography scans (CT), magnetic resonance images (MRI) and MR diffusion tensor images (MR-DTI). Since the resolution of these observations is limited, typically around 1mm x 1mm x 1mm in the best case, observable factors are limited. Mathematical models at the macro-scale try to formulate the tumour growth using observations coming from this scale. For this reason, these models include fewer factors and are mathematically simpler than the microscopic models of Section 1.2.1. On the other hand, while microscopic models simulate the tumour growth in theoretical settings (infinite boundaries, known location of different structures,...), macroscopic models use real settings personalised for each patient , e.g. real boundaries of the brain, grey-white matter segmentation...

Based on the effect of the tumour on the brain, macroscopic models can be classified in two different classes: mechanical models, which concentrate on the mass-effect of the tumour on the brain tissue, and diffusive models, which concentrates on the invasion of surrounding tissue by tumour cells. From the mathematical point of view, all macroscopic models use continuum formulations based on a continuous local density or proportion of tumour cells. As a result, formulations contain several ordinary and/or partial differential equations to describe the growth process.

• Diffusive Models:

Almost all diffusive macroscopic models use the reaction-diffusion formalism. This formalism models the invasive tumour by adding a diffusion term to simple solid tumour growth models, which formulates proliferation of cells, see table 1. The "building block" equation of this formalism is the partial differential equation (PDE) given as:



where u is the tumour cell density, D is the diffusion tensor of tumour cells and R(u,t) is the so-called reaction term. This equation isolates two different characteristics of the tumour growth in the two terms: diffusion and proliferation. The first term on the right hand side describes the invasion of tumour cells by means of a directed Brownian motion, which is characterised by the diffusion tensor D. The second term in the equation, R(u,t), describes the proliferation of tumour cells. The second equation represents the no-flux boundary condition which is applied at the skull and ventricles, formulating the fact that tumour cells do not diffuse towards these structures.

Recent works on diffusive models concentrates on the infiltration part of the equation, in other words on the diffusion tensor. Following the experimental results showing the high infiltration speed of tumour cells on the myelin sheath, several models proposed to use an isotropic but heterogeneous D varying depending on the tissue (grey-white). Latest works also included the hypothesis that tumour cells actually follow the fibre tracts by creating an anisotropic tensor D.

• Mechanical Models:

Mechanical models, which concentrated on the mass-effect of the tumour, contain two distinct coupled formulations, one for the tumour growth and one for the mechanical characteristics of the brain tissue. These model tumour growth using rather a simpler formulation consisting only of mitosis (population growth model) causing an outwards pressure. They concentrate on the mechanical interaction between a growing tumour and the brain tissue whose characteristics are modelled with detail. There have been many works on characterising the mechanical properties of the brain tissue, which is deformable but not elastic.

In [Wasserman and Acharya, 1996] it is said that the brain tissue is a sponge like material, possessing instantaneous properties of elastic materials and time-dependent properties of the viscoelastic ones. Moreover, there is a great variation between elastic parameters of brain tissue within similar tissues as well as between differing tissues. Instead of formulating these complex mechanical characteristics, almost all models use assumptions to simplify brain tissue's characteristics. First works used linear elastic material properties for the brain tissue simplifying the model. However, this assumption did not allow large deformations usually observed in clinical cases. More recent works relaxed this assumption and used hyper elastic material properties obtaining these large deformations.

On the other hand, some recent works kept the linear elastic material assumption but tried to improve the tumour growth model. Eventually they obtained a coupling between more complicated diffusive models explained above and a simpler mechanical model describing the deformation of the brain tissue.

1.3 References

• [Araujo and McElwain, 2004] Araujo, R. and McElwain, D., 2004. A history of the study of solid tumor growth: the contribution of mathematical modeling. Bul. Math. Biology, 66.
• [Mantzaris et al., 2004] Mantzaris, N., Webb, S., and Othmer, H., 2004. Mathematical modeling of tumor-induced angiogenesis. Journal of Math. Biology, 49.
• [Orme and Chaplain, 1996] Orme, M. and Chaplain M.A.J., 1996. A mathematical model of vascular tumor growth and invasion. Math. Comp. Modelling, 23.
• [Wasserman and Acharya, 1996] Wasserman, R. and Acharya, R., 1996. A patient-specific in-vivo tumor model. Math. Biosci.,136.
• Zheng, X., Wise, S., and Cristini, V., 2005. Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method. Bull. Math. Biology, 67.

Tumour Growth Modelling in the Brain

2. A Macroscopic Patient Specific Growth Model

The aim of mathematical growth models is to explain the underlying dynamics of the tumour growth in detail. As a result these models can be very complicated and very specific to certain situations. Patient specific modelling on the other hand, aims at creating flexible models that can be adapted to different patient cases, with obviously a trade-off in the level of detail included. This ability to adapt is crucial to build the link between theoretical models and clinical applications. Such models are usually less complicated and have few parameters that can be measured from clinical data. The biggest advantage of patient specific models is that they can be used in therapy planning and prediction of the growth for the specific cases.

In patient specific growth modelling, the formulations are based only on the observations that can be gathered in the clinical environment such as images, genetic data and histology. Here, we present a patient specific model based on MR images proposed in [Clatz et al., 2005]. This model concentrates on the information coming from the anatomical (T1 and T2 weighted) images and diffusion weighted images. T1 and T2 images are used to segment/delineate the bulk tumour (gross tumour volume 1 - GTV1), the tumour infiltrated oedema (gross tumour volume 2 - GTV2) and the white matter. Diffusion weighted images are used to extract fibre directions, which guides the tumour infiltration (See the figure below). The growth model is divided into two different parts: the diffusion, and the mechanical model.

Data used in the model and the information extracted: (a) T1 and T2 weighted images of a patient suffering from a high grade glioma. (b) The diffusion tensor images derived from the diffusion weighted images. (c) Different tissue structures extracted from the anatomical images like, the white matter (3), the grey matter (2), the falx cerebri (5), the skull (1) and the ventricles (4) are embedded in a mesh. (d) The bulky part of the tumour consisting the necrotic core (1) (if it exists). (e) The tumour infiltrated oedema region (2).

2.1 Diffusion Model

The diffusion model explains the infiltration of tumour cells into the brain parenchyma and the proliferation of tumour cells. The model is formulated in the reaction-diffusion formalism as given below. Change of tumour cell density at a point u(x) is described by an anisotropic diffusion and a nonlinear reaction process.

The infiltration of tumour cells is explained by the diffusion process ∇(D(x)∇u), which is characterised by the diffusion tensor D. The proliferation of tumour cells are embedded in the reaction part ρu(1-u) with the mitosis rate ρ given in terms of doubling rate. The Neumann boundary conditions (D∇u)n dictates that tumour cells will not pass through the skull nor the ventricles. The diffusion tensor D is constructed following the observation of [Giese et al. 1996] saying that tumour cells move faster on the myelin sheath and also discussions with neurosurgeons who suggested that tumour cells follow the fibre tracts. Thus, D is an isotropic tensor in the grey matter infiltrating in all directions equally. In the white matter, D is an anisotropic tensor following the white fibre tract with a higher mean diffusion than in the grey matter:

The following figure summarises the diffusion properties of different tissues in the brain.

Diffusion properties of different tissues in the brain. Tumour cells diffuse through the fibre tracts while they do normal diffusion in the grey matter with a lower speed. The boundary conditions dictate the cells cannot pass through the brain boundary.

2.2 Mechanical Model

The mechanical model concentrates on the mass effect of the tumour growth on the surrounding structures. As the tumour grows, due to abnormal speed of mitosis and the infiltration, surrounding tissues are deformed. This deformation is formulated by coupling the pressure caused by the change in tumour density with the material characteristics of the brain structures. Brain parenchyma is an extremely complicated material and a full mathematical description would be complex. Here, the parenchyma is assumed to be a linearly elastic material for simplicity. The constitutive equations for the brain parenchyma is given with the following system.

where ε is the stress tensor defining the pressure exerted by the tumour growth, σ is the strain tensor defining the deformation, λ and μ are the Lamé coefficients and v is the displacement. Other material properties are summarised in the following figure.

Material properties of different structures in the brain. The parenchyma is assumed to be a linearly elastic material. The ventricles are modelled as a free body. The falx cerebri is also a linearly elastic material but much more stiffer than the parenchyma. Lastly the skull is modelled as a non-elastic structure.

The pressure exerted by the growth of the tumour is modelled in two different parts. The bulky part of the tumour grows solely by proliferation without diffusion. The mechanical effect of this proliferation is modelled as a constant outward pressure describing expansion of GTV1. On the other hand, the mass effect of the infiltrating tumour cells is described by a local pressure based on the tumour density. The mass effect formulation is given as follows.

where F_e represent the external forces and the α is the coupling factor between the diffusion model and the mechanical model.

2.3 Growth Simulation

Using the models explained above one can simulate the tumour growth for a specific patient case. The output of the model for a patient suffering from a high grade glioma is shown in the following video.

Using two different images of the patient taken at different times we can compare the simulated growth with the actual growth. We observe that the simulated growth highly resembles the actual observed growth between the successive MR images although there is still a difference. This difference can be reduced by optimising the parameters of the model for the patient and using better segmentation techniques.

First and the third columns show MR images taken four months apart from a patient with high grade glioma. Second column show the initial condition for the model based on the initial delineation of the tumour. Fourth column shows iso-contours of tumour cell density computed by the growth model.

2.4. References

• [Giese et al., 1996] Giese, A., Kluwe, L., Laube, B., Meissner, H., Berens, M., Westphal, M., 1996. Migration of human glioma cells on myelin. Neurosurgery, 38.
• [Clatz et al., 2005] Clatz, O., Sermesant, M., Bondiau, P. , Delingette, H., Warfield, S. Malandain, G., Ayache, N., 2005. Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation. IEEE TMI, 24.

Tumor Growth Modelling in the Brain

3. An Eikonal Approximation for the Reaction-Diffusion Type Tumor Growth Models

The reaction-diffusion type tumor growth models as explained here are able to describe the macroscopic growth of gliomas as visible in medical images. Although these models are successful in simulating growth of realistic tumors there is one important discrepancy between these models and the medical images. Reaction-diffusion type growth models formulate the evolution of tumor cell density distribution at every point in the brain. However, tumor cell density distribution is not available in clinically used medical images. We observe delineations around the tumor or around the oedema depending on the modality of the image. In accordance with this, the growth of the tumor observed from the images corresponds to the evolution of these delineations rather than the evolution of tumor cell densities.

MR Flair images of a second grade astrocytoma:
(Left) image at the first examination (Right) image at the second examination.
In the anatomical MR images we observe the boundary
of the visible part of the tumor (the tumor delineation) rather than the tumor cell densities.

This discrepancy becomes crucial especially when one wants to extract model parameters from the images, in other words when solving the inverse problem. One solution to this problem might be to assume tumor cell density distributions based on the visible delineation. However, this adds extra assumptions about the image acquisition process. Our approach is to avoid this extra assumption and derive a reduced model. The definition of the reduced model is that it will formulate the evolution of the tumor delineation in time based on the same growth dynamics as the reaction-diffusion model. However, it will not use the information of tumor cell density distribution while doing so. Using the asymptotic properties of the reaction-diffusion equation we derive such a reduced model.

3.1. Asymptotic Properties of Reaction-Diffusion Equations

Get the Flash Player to see this player.

Convergence of Reaction-Diffusion Evolution in the Infinite Cylinder: The PDE is initialized with a compact support function and let evolve in time on the left. It reaches its asymptotic shape, the traveling wave. On the right we follow the speed of the u=0.5 point, the blue dot on the left. Observe that the speed of this point approaches the asymptotic speed.

The speed of the traveling wave (shown in dark)
is time dependent and it converges to the asymptotic
speed (shown in dashed). The convergence is slow
therefore, we need to take it into account if we want
to approximate the speed of this wave.
The approximation we use takes this into account
(shown in dashed-dots).

3.2. Eikonal Approximation for Tumor Growth

Comparison of the motion of the tumor delineation computed through simulations with reaction-diffusion model and thresholding the cell densities (in white) and the reduced model without using tumor cell densities (in black). We observe that the high similarity between the images. Each contour represents the position of the tumor delineation at a single time point from day 400 to 1200 as we move outwards.

3.3. References

• [Konukoglu et al., 2007a] Konukoglu, E., Sermesant, M., Clatz, O., Peyrat, JM, Delingette, H., Ayache, N., 2007. A recursive anisotropic fast marching approach to reaction diffusion equation: application to tumor growth modeling., IPMI 2007.
• [Ebert and Saarlos, 2000] U. Ebert, W.V. Saarloos, 2000. Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts. Physica D: Nonlinear Phenomena, 146.
• [Aronson and Weinberger, 1978] ]D.G. Aronson and H.F. Weinberger, 1978. Multidimensional Nonlinear Diffusion Arising in Population Genetics. Advances in Mathematics, 1978, 30.

Tumour Growth Modelling in the Brain

4. Extrapolating Tumour Invasion

The diffusive nature of glial tumours inhibits the apparent boundary between cancerous and healthy regions, decreasing the effects of conventional therapies like surgical resection and radiotherapy. The current imaging techniques are able to detect only a part of the total tumour [Tracqui et al., 1995]. In the literature, there are several suggestions for the lowest detection limit of CT images in terms of number of tumour cells per unit volume or area. The Figure following this paragraph illustrates the tumour profile which is visible in the image (the solid line curve), and the theoretical total tumour profile. The usual clinical practice is to remove surgically the visible part, and to try to kill tumour cells which are not visible in the images using radiotherapy and/or chemotherapy. The conventional radiotherapy treatment uses a geometric approach and applies radiation in a volume defined by constant margin of about 2 cm around the surgical site. This approach does not take into account the different speeds of tumour diffusion for the grey and the white matter [Giese, 1996]. Using the same margin for both tissues may kill healthy grey tissue while not reaching the extent of tumour cells in the white matter, leaving alive cancerous cells on the fibres. By constructing physiologically determined radiotherapy margins taking into account the tumour growth dynamics and the properties of the tissue, the therapy may be more effective.

Tumour cell distribution is assumed to be continuous and smooth as opposed to the signal intensity we see in medical images.

The work explained here proposes a method to estimate the distribution of tumour cells for low density regions using a single image and a set of estimated parameters [Konukoglu et al, 2006]. Our formulation constructs an approximation of the tails of the tumour distribution (invasion margins) which are not visible in images. In addition, this formulation creates appropriate initial conditions for tumour growth models, which can be used for simulating tumour growth and predicting the distribution of the tumour in a future time. Using an initial condition with well estimated tails will increase the ability of these models to predict not only the tails of the tumour density, but also the whole distribution.

4.1 Method

Our formulation is based on the Fisher-Kolmogorov equation, which is been widely used to model the growth of brain tumours, [Murray, 2002]:

where u is the tumour cell density, ρ is the average proliferation rate of tumour cells and D is the diffusion tensor which characterises the diffusion behaviour of cancerous cells. D is taken to be spatially varying with a higher determinant on the white matter yielding a higher diffusion rate. It is also taken to be an anisotropic tensor following the fibres in the brain using the diffusion tensor images (DTI), as proposed in [Clatz et al., 2005]. Using this setting we derive the formulation describing the profile of the tumour front at low cell density parts at a single time instant:

where u₀ is the lowest cell density threshold for images and Γ is the last visible contour of the tumour in the image. This derivation is based on the fact that Fisher-Kolmogorov equation admits a travelling wave solution in 1D. Using this formulation, given a single image with a tumour and a set of estimated parameters (D and ρ), we are able to extrapolate the unseen parts of the tumour.

4.2 Invasion Margins

The advantage of this formulation is that it extrapolates the low density profile of the tumour using a single time instant, providing a tool for medical doctors to use in planning radiotherapy. Moreover, the formulation constructs a patient-specific tumour profile based on the underlying tissue structures and the current location of the tumour, which provides a realistic initial condition for tumour growth models in trying to predict the evolution of the tumour.

Left: A virtually grown tumour segmented in red and plotted with the 2cm constant radiotherapy margin around it. Middle: Extrapolated distribution of tumour cells beyond the visible part of the tumour. Right: Comparison of tumour cell distribution and constant therapy margin. Red part is the successful region of irradiation while brown and green parts are not. We observe that the tumour invasion might be different than the conventional margin.

In the Fisher-Kolmogorov formulation of tumour growth, speed of growth depends mostly on the product of D and ρ, while the extent of the tumour depends on the ratio of these parameters. When we have two images of the same patient taken at different time instances, the observed evolution of the image tumour delineation may be explained by an infinite number of parameter since the same product will produce very similar evolutions of the front wave. However, the extent of the infiltration margin will be very different as it depends on the ratio of the parameters. This effect is shown in the video below.

4.3 References

• [Clatz et al., 2005] Clatz, O., Sermesant, M., Bondiau, P., Delingette, H., Warfield, S. Malandain, G., Ayache, N., 2005. Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation. IEEE TMI, 24.
• [Giese et al., 1996] Giese, A., Kluwe, L., Laube, B., Meissner, H., Berens, M., Westphal, M., 1996. Migration of human glioma cells on myelin. Neurosurgery, 38.
• [Konukoglu et al, 2006] Konukoglu, E., Clatz, O., Bondiau, P., Delingette, H., Ayache, N., 2006. Extrapolating tumor invasion margins for physiologically determined radiotherapy regions. Proc of MICCAI'06, LNCS 4190, 2006.
• [Murray, 2002] Murray, J., 2002. Mathematical Biology. Springer-Verlag.
• [Tracqui et al., 1995] Tracqui, P., Cruywagen, G., Woodward, D., Bartoo, G., Murray J., Alvord, E., 1995. A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth. Cell Prol. 28.
• Swanson, K., Alvord, E., Murray, J., 2002. Virtual brain tumors (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy. Br. J. Cancer, 86.
• Jbabdi, S., Mandonnet, E., Duffau, H., Capelle, L., Swanson, K., Pelegrini-Issac, M., Guillevin, R., Benali, H., 2005. Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging. Mag. Reson. in Med., 54.

Tumour Growth Modelling in the Brain

5. Identifying Model Parameters

In cancer treatment, understanding the aggressiveness of the tumour is essential in therapy planning and patient follow-up. In this work we propose a method for quantifying the progression of the critical target volume (CTV) of glial-based tumours, on the basis of their growth dynamics. The formulation is based on the tumour growth model proposed by [Clatz et al., 2005], which uses reaction-diffusion formalism:

With the proposed method, we obtain quantitative estimates for the speed of invasion in white and grey matter by solving the patient specific parameter identification problem for this growth model using MR images taken at two different time instances, t₁ and t₂, from the same patient. The parameter identification problem is formulated using the front approximation of reaction-diffusion equations, which results in anisotropic Eikonal equations. The anisotropic fast marching method proposed in [Konukoglu et al., 2007a] is used for numerical solutions yielding an efficient algorithm.

5.1 Method

The model given above requires tumour cell density u to be known at every point as an initial condition. However, this is not the case for medical images where only contours around gross tumour volume (GTV) and CTV are available. The front motion approximation of reaction-diffusion equations offers a solution for this discrepancy between information needed and observations available [Konukoglu et al., 2007b]. Taking the contour around the CTV, Γ, as the last visible tumour front we can use such an approximation to model its evolution. As a result we obtain the following travelling wave formulation:

where T(x) represents the time at which Γ passes from point x.

In order to formulate the inverse problem we follow the modelling assumption as given in [Clatz et al., 2005] and state T(Γ₁)=0 and T(Γ₂)=t₂-t₁, where Γ₁ and Γ₂ corresponds to contours around CTV regions observed in images taken at time t₁ and t₂ respectively. Parameter identification process tries to find parameters d_g and d_w that create a T function that would satisfy these conditions. Notice that ρ is a multiplicative factor in the travelling time formulation and it cannot be determined independently from the D matrix by just looking at the motion of the tumour front. To tackle this, we treat ρ as a known constant in the parameter identification problem. t₁ is not available in clinical circumstances hence we use t₂-t₁. The formulation can be given as the minimisation problem

where C is the objective function to minimise with respect to d_w and d_g, Γ₂^c is the computed contour using the front approximation with the given parameters and dist(A,B) is the distance between two iso-surfaces taken as mean distance from voxels of A to the closest voxel of B.

Left: T2 weighted MR images of a patient taken at t₁ and t₂ and the tumour delineations around the tumour infiltrated oedema Γ₁ and Γ₂ manually done. Right: The minimisation surface C(d_w, d_g) for the parameter identification problem given by the images on the left.

5.2 Minimising Parameters

The method explained in the previous section allows us to identify the diffusion related parameters of the reaction-diffusion based growth model. In identifying these parameters we solve the quantification problem: How fast is the tumour growing? The way we answer this question is by providing two parameters corresponding to speed of diffusion of tumour cells in the grey matter and in the white matter. Such an approach quantifies the spiky nature of invasion of brain tumours, which can be very hard to capture using methods such as measuring the diameter of the tumour. Although we demonstrate the method on high grade gliomas, same idea extends to low grade gliomas which are modelled using similar formulations.

Sagittal and axial slices of a patient suffering from a high grade glioma is given. The underlying image is taken at t₂. Blue contours: The tumour delineation at t₁ manually done by an expert, Γ₁. Black contours: The tumour front computed by the travelling time formulation Γ₂* with the parameters minimising the functional d_w* and d_g*. The minimising speed values in the grey and in the white matter are also given in the image. Observe how the computed black contour is close to the CTV in the underlying image. Observe the very high difference between white matter and grey matter speed which is coherent with the experimental observations.

5.3 References

• [Clatz et al., 2005] Clatz, O., Sermesant, M., Bondiau, P., Delingette, H., Warfield, S. Malandain, G., Ayache, N., 2005. Realistic simulation of the 3d growth of brain tumors in mr images coupling diffusion with biomechanical deformation. IEEE TMI, 24.
• [Konukoglu et al., 2007a] Konukoglu, E., Sermesant, M., Clatz, O., Peyrat, JM, Delingette, H., Ayache, N., 2007. A recursive anisotropic fast marching approach to reaction diffusion equation: application to tumor growth modeling., IPMI 2007.
• [Konukoglu et al., 2007b] Konukoglu, E., Clatz, O., Bondiau, P., Sermesant, M., Delingette, H., Ayache, N., 2007. Towards an Identification of tumor growth parameters from time series of images, to be presented at Miccai 2007.

Tumor Growth Modelling in the Brain

6. Monitoring Slowly Evolving Tumors

Most of the pediatric brain tumors as well as the majority of adult tumors are slowly evolving pathologies. Although modeling the growth of these tumors is important the slow rate of evolution brings up another question: detecting the growth. More explicitly, the clinically more relevant question is detecting the change of tumor volume between two images of a patient taken at different times. Here we describe an approach thatsemi-automatically performs this task using longitudinal medicalimages. We specifically focus on meningiomas, which expertsoften find difficult to monitor as the tumor evolution can beobscured by image artifacts.

Two scans of a pediatric patient with a diffusive glioma.
The scan on the left is taken several months before the scan on the right.
The evolution of the tumor is slow therefore detecting it becomes
a crucial question.

Meningiomas are the most common type of primary brain tumor. Most ofthese tumors are categorized as benign pathology that grows slowlybetween brain tissue and dura. To avoid the risk of surgery,neurosurgeons carefully monitor patients with benign meningiomas byhaving the patient regularly undergo Magnetic Resonance (MR)scanning. An expert assess the tumor growth through visual inspectionof consecutive 3D scans. A precise analysis, however, isextremely difficult as slow growth is often obscured bychanges in the head position or intensity profile between the twoscans. We address this issue by describing a relatively fast androbust method that semi-automatically analyzes the tumor evolution.

Two scans of a meningioma patient. The scan on the left
is taken several months before the scan on the right. The pose and intensity
differences makes it hard to detect small volume changes.

Our approach first semi-automaticallysegments the tumor in the initial patient scan. It then aligns thesecond scan of the patient to the first using a hierarchical rigidregistration approach. Finally, it measures growth or shrinkage fromthese images, for which we suggest two different types ofmetric. Motivated by [Angelini et al. 2007], the first metric detects change through analyzingdifferences in intensity distributions. We differ from [Angelini et al. 2007] in that werelate the analysis to hypothesis testing. The second metricis motivated by the work of [Rey et al. 2002], which detects change byanalyzing the deformation field between two scans. However,our metric returns a quantitative analysis of thedifferences in volume (mm³). Our approach also differssignificantly from works on segmentation, suchas [Prastawa et al. 2003, Liu et al. 2005] as we estimate volume change by simultaneously analyzing the sequence of scans.

6.1. The Framework

The software pipeline is composed of three steps, which are:

• tumor segmentation,
• image registration and
• change detection.

First step is to perform the semi-automatic segmentation
of the first scan. User selects a region of interest and
determines the right intensity. The algorithm applies
thresholding and island removal to obtain the final
segmetation.

The semi-automatic segmentation is based on a user-defined boundingbox around the tumor and a lower bound of the intensitiesthat characterize the pathology. From these indicators, the pipeline canreliably extract most of the pathology when it appears as homogeneous, bright objects in the MR images. Most of the slowly growing tumors such as low grade astrocytomas and meningiomas appear as such objects. For more complicated pathologies different segmentation techniques may be utilized and these can be plugged into the pipeline easily. The pipeline post-processes the resulting binary map by removing small islands and holes caused by the noise in the MRs. We note that the resulting map mayalso include part of the dura and vessels since these structures may haveintensity patterns that are similar to pathology. However, becausethey are static, these additional structures should not substantially impact the analysis.

The two scans may be misaligned due to different head pose
during the image acquistion. The second step in our pipeline
is to rigidly align the two scans. The image on the left
shows the bounding box and the segmentation we have obtained
from the first scan drawn on to the second scan.
Notice the misalignment. The image on the left shows
the same thing after the rigid alignment.

The second step automatically aligns the pathology ofthe second scan to the first. It does so by rigidly registering the scans. Now, we assume that the previous bounding box islarge enough so that its coordinates also define a bounding box around thepathology in the second scan. We then accommodate partial voluming inimages by increasing the resolution of both bounding boxes. Finally, theframework addresses non-linear perturbation artifacts caused by the MRacquisition by rigidly aligning the contents of the bounding boxes with each other. This results in two images where, in theory, barring temporal changes, the pathology is well aligned.

The final step of our approach measures the tumor evolution based onthe initial segmentation and the bounding boxes described above. We now proposetwo different types of analysis for detecting the volume change.

After the registration step we obtain two well-aligned
images of the pathology in different times, in the
upper row we show these images. At this point the
images are ready to be compared for growth detection.
The image on the bottom shows a map of volume
changes obtained using the Intensity Analysis.

• Intensity Analysis: The first analysis is based on a statistical model that differentiates changes in the intensity patterns that are due to evolving pathology rather than image artifacts. In this model we first learn the intensity variability between two scans that is due to image artifacts, σ_art. This is simply achieved by comparing the intensity patterns in the two scans using the regions where there are no enhancement due to the tumor. Based on the learnt σ_art we then compare the intensity values in the scans and mark the differences exceeding σ_art as growth or shrinkage.
• Deformation Analysis: The second analysis is based on the amount of deformation required to map the bounding box coming from the first scan to the one coming from the second scan. We obtain this deformation field, D, by non-rigidly registering the two scans using the diffeomorphic demons algorithm [Vercauteren et al. 2007]. Using the D we can estimate the amount of volume change by two ways. The first one is to analyze the Jacobian of this field J_D whose determinant provides us the volume change. The second way is to apply the deformation to the delineation of the first scan obtained using semi-automatic segmentation. Comparing the original and the deformed delineations gives the us the volume difference.

6.2. The Output

	Intensity Analysis	Jacobian Analysis	Segmentation Analysis
Computed Growth	557.4 mm3 (7.8%)	1539.4 mm3 (21.6%)	1555.7 mm3 (21.9%)

The two unaligned and unnormalized scans seen above are compared in a robust manner
using the framework explained in the previous section. The volume differences
are given in mm³ and in %.

6.3. References

• [Konukoglu et al. 2008] E. Konukoglu, W. M. Wells, S. Novellas,N. Ayache, R. Kikinis, P. M. Black, K. M. Pohl. 2008, Monitoring Slowly Evolving Tumors. ISBI, pp. 812-815.
• [Angelini et al. 2007] E.D. Angelini and J. Atif and J. Delon and E. Mandonnet and H. Duffau and L. Capelle. 2007, Detection of glioma evolution on longitudinal {MRI} studies. ISBI, pp. 49-52.
• [Rey et al. 2002] D. Rey and G. Subsol and H. Delingette and N. Ayache, 2002, Automatic detection and segmentation of evolving processes in 3D medical images: Application to multiple sclerosis. Medical Image Analysis, vol. 6, pp. 163-179.
• [Prastawa et al. 2003] M. Prastawa and E. Bullitt and N. Moon and K.V. Leemput and G. Gerig, 2003, Automatic brain tumor segmentation by subject specific modification of atlas priors. Academic Radiology, vol. 10.
• [Liu et al. 2005] J. Liu and J. Udupa and D. Odhner and D. Hackney and G. Moonis, 2005, A system for brain tumor volume estimation via {MR} imaging and fuzzy connectedness. Computerized Medical Imaging and Graphics, vol 29.
• [Vercauteren et al. 2007] T. Vercauteren and X. Pennec and A. Perchant and N. Ayache, 2007, Non-parametric diffeomorphic image registration with the demons algorithm. MICCAI, vol. 4792, pp. 319-326.

Computational Anatomy of the Brain

1. Computational Anatomy: Aims and Methods

1.1 Evolution and Revolutions in Anatomy

According to the dictionary, anatomy is the science that studies the structure and the relationship in space of different organs and tissues in living systems. During the antiquity and middle-ages, anatomical descriptions were mainly based on animal models and the physiology was more philosophical than scientific. Modern anatomy really began at the renaissance with the authorisation of the dissection of human cadavers, giving birth to the "De humani corporis fabrica" published by in 1543 by Vesale (1514-1564), and was strongly pushed by the progresses in surgery, as exemplified by the "Universal anatomy of the human body" (1561-62) of the great surgeon Ambroise Paré(1509-1590).

Figure 1. Left: first cerebral atlas, Vesale, 1543. Middle: human anatomy (skeleton) by Paré, 1585. Right: Atlas of the brain by Talairach and Tournoux.

During the following centuries, many progresses were done in anatomy thanks to the introduction of new observation tools like microscopy and histology, going down to the level of cells in the 19th and 20th centuries. However, in-vivo and in-situ imaging is radically renewing the field since the 1980ies. An ever growing number of imaging modalities allows to observe both the anatomy and the function at many spatial scales (from cells to the whole body) and at multiple time scales: milliseconds (e.g. beating heart), years (growth or ageing), or even ages (evolution of species). Moreover, the non-invasive aspect allows repeating the observations on multiple subjects, giving rise to large databases. This has a strong impact on the goals of the anatomy which are changing from the description of a representative individual to the description of the structure and organisation of organs at the population level. The huge amount of information generated (typically 50 to 150 Mb for a clinical MRI exam, much more for a CT one) also raises the need for computerised methods to extract and structure information. This led in the last 10 to 20 years to the gradual evolution of descriptive atlases into interactive and generative models, allowing the simulation of new observations. Typical examples are given for the brain by the MNI 305 ([Evans et al, 1993]) and ICBM 152 ([Mazziotta et al, 2001] ) templates that are the basis of the Brain Web MRI simulation engine [Collins et al, 1998]. In the orthopaedic domain, one may cite the "bone morphing" method [Fleute and Lavallée, 1998,Rajamani et al, 2004] that allows to simulate the shape of bones.

1.2 Aims of Computational Anatomy

The combination of these new observation means and of the computerised methods is at the heart of computational anatomy, an emerging discipline at the interface of geometry, statistics and image analysis whose goal is to develop algorithms to model and analyse the biological shape of tissues and organs. The goal is not only to estimate representative organ anatomies across species, populations, diseases, ageing, ages but also to model the organ development across time (growth) and to establish their variability. Another goal is to correlate this variability information with other functional, genetic or structural information (e.g. fibre bundles extracted from diffusion tensor images). From an applicative point of view, a first objective is to understand and to model how life is functioning at the population level, for instance by classifying pathologies from structural deviations (taxonomy) and by integrating individual measures at the population level to relate anatomy and function. For instance, the goal of spatial normalisation of subjects in neuroscience is to map all the anatomies into a common reference system. A second application objective is to provide better quantitative and objective measures to detect, understand and correct dysfunctions at the individual level in order to help therapy planning (before), control (during) and follow-up (after). This generally implies mapping some generic (atlas-based) knowledge to patients-specific data through atlas-to-patient registration. Computational anatomy is currently a very active research field, as exemplified by the IPAM's Graduate Summer School on Mathematics in Brain Imaging organised by P. Thompson, M. Miller in 2004 or the workshop on the Mathematical Foundations of Computational Anatomy (MFCA'06), organised in September 2006 in conjunction with MICCAI'06, the 9th International Conference on Medical Image Computing and Computer Assisted Intervention.

1.3 Methods of Computational Anatomy

In the case of observations of the same subject, many geometrical and physically based registration methods were proposed to faithfully model and recover the deformations. However, in the case of different subjects, the absence of physical models relating the anatomies leads to a reliance on statistics to learn the geometrical relationship from many observations. The general method is to identify anatomically representative geometric features (points, tensors, curves, surfaces, volume transformations), and to model their statistical distribution across the population. This can be done, for instance, via a mean shape and covariance structure after a group-wise matching. In the case of the brain, one can rely on a hierarchy of structural models ranging from anatomical or functional landmarks like the AC and PC points ([Talairach and Tournoux, 1988; Bookstein, 1978]), curves (e.g. crest lines, sulcal lines) ([Mangin et al 2004; Le Goualher et al. 1999, Fillard et al. 2005]), surfaces like sulcal ribbons ([Thompson et al., 1997; Andrade et al., 2001; Vaillant et al., 2007]), images seen as 3D functions, which lead to voxel-based morphometry (VBM [Ashburner and Friston, 2000]), diffusion imaging or rigid, multi-affine or diffeomorphic transformations ([Trouvé 1998; Miller et al., 2002; Arsigny et al., 2006]), leading to Tensor-based morphometry (TBM).

Figure 2. Top: Surface of the cortex extracted from 4 different subjects with the Anatomist/BrainVISA package. Bottom: example extraction and labelling of the sulci with the same tools. One can see that finding the invariants of the cortex shape is a difficult problem. (Courtesy of J.-Fr. Mangin, CEA, Laboratoire de neuroimagerie (LNAO), Neurospin, Saclay, France)

One crucial point is that these features usually belong to curved manifolds rather than to Euclidean spaces, which precludes the use of classical linear statistics. For instance, the average of points on a sphere is located inside the sphere and not on its surface. To address this problem, one has to rely on statistical tools that work directly on Riemannian manifolds in an intrinsic way. A consistent set of such tools including the mean value, covariance matrix, Normal law, Mahalanobis distance were for instance proposed in [Pennec, 2006] along with efficient algorithms and tractable approximations for small variances. These tools were later extended to a complete algorithmic framework to compute on manifolds with the example of so-called "tensors" (symmetric positive definite matrices) [Pennec et al, 2006]. In particular, it was showed that one can perform interpolation, filtering, isotropic and anisotropic regularisation and restoration of missing data (extrapolation or in-painting) on manifold valued images by using generalised weighted means and partial derivative equations (PDEs).

1.4 Example use of Computational Anatomic Models in the Clinical Workflow

1.4.1 Better Constrain Atlas-to-Patient Registration for Radiotherapy

The planning step for conformal radiotherapy requires the accurate localisation of the tumour, to maximise its irradiation, and of the critical structures where the irradiation has to be minimises. To segment these structures for each patient, one standard method is to register a previously labelled atlas to the patient image. This allows transferring the generic atlas segmentation toward the patient-specific space. This segmentation can then be used directly, or as an initialisation for a more complex segmentation algorithm [Commowick, 2007]. In such a system, the main difficulty is to obtain an inter-subject registration algorithm which is accurate enough and, more importantly, robust to the anatomical variability and to the pathologies (tumours may be quite important).

Figure 3. Structures of interest for conformal radiotherapy segmented by P.Y. Bondiau at CAL, Nice based on the MNI brain atlas [Bondiau, 2004].

The main method to enforce meaningful deformations is to penalise non sensible ones through a regularisation criterion. Some authors used physical models like elasticity or fluid models [Bajcsy and Kovacic, 1989, Christensen et al., 1997]. For efficiency reasons, other authors proposed to use non-physical but fast regularisation methods like Gaussian filtering [Thirion, 1998, Pennec et al., 1999, Modersitzki, 2004]. However, since we do not have in general a model of the deformation of organs across subjects, no regularisation criterion is obviously more justified than the others. Thus, most of the existing work in the literature rather tries to capture the organ variability from a statistical point of view on a representative population of subjects (see e.g. [Thompson et al., 2000, Rueckert et al., 2003, Fillard et al., 2005]). For instance, following [Lester et al., 1999], the algorithm RUNA developed in [Stefanescu et al., 2004, Stefanescu, 2005] was a first attempt to obtain a computationally efficient but highly steerable nonlinear registration algorithm that includes some anatomical information about the tissue types. It uses a non-stationary transformation regularisation which is strong where the local deformability is expected to be low, and conversely. This also allows taking into account pathologies such as tumours or previous resections.

Figure 4. Registration of the atlas towards the patient image to segment structures at risk for radiotherapy planning. a) Patient image, b) Stiffness field, c) Atlas deformed towards the patient images without taking into account the resection, d) Atlas deformed towards the patient images with the pathology taking into account.

Moreover, the regularisation can be locally tuned along spatial directions through the use of a tensor field, as we did for instance in [Commowick et al., 2005], by introducing a method to compute scalar and tensor based deformability statistics over a database of patient.

Figure 5. Comparative results of the atlas-based segmentation using RUNA guided by the heuristic scalar stiffness map (left), the scalar statistical stiffness map (middle) and the statistical tensor stiffness map (right). Top: image of the stiffness map used. Bottom: Sagittal slice of the 3D patient image with the segmentation superimposed. The inclusion of the statistical information clearly allows to better deform the brain stem area in an anatomically more meaningful way, resulting in a better and smoother segmentation. Brain images are courtesy of Dr. P.-Y. Bondiau, CAL, Nice.

1.4.2 Mapping the fibres to patient images for tumour growth models

As shown in the tumour growth model section, the DTI information is crucial in the formulation of the tumour growth model because tumour cells move much faster on the white matter and they tend to follow the fibre tracts. Hence, the infiltration pattern of the tumour highly depends on the local fibre directions. Since diffusion tensor images give the average local fibre orientation they can be used to extract these tumour highways as it was done in the model we have explained to construct the diffusion tensor of tumour cells. Unfortunately, DT-MRI images are seldom acquired in the clinical routine because of technical constraints (old clinical MR scanners often do not have gradients powerful enough to acquire this sequence), time constraints which are crucial for the image quality, especially on children (a typical DTI acquisition lasts about 10 minutes during which the child should not move, which may require sedation), ethical issues and economical efficiency (there is no clinical indication for DTI yet).

Figure 6. Example mapping of tensors from a healthy subject to a patient. (a) T1 weighted image of the healthy subject (atlas). (b) T1 weighted image of the patient. (c) DTI of the healthy subject in the standard colour mode. The colour encodes the principal direction of the tensor (green: front/back, red: left/right, blue: top/down), while the intensity is related to the fractional anisotropy which quantifies the reliability of this principal direction. (d) Result of the resampling of image (c) using the non-linear transformation found by the registration of the anatomical images (a) and (b).

In such a personalisation process, the quality of the registration algorithm is essential as it directly affects the model and the growth simulation. Thus, using prior knowledge about the inter-subject variability is necessary to obtain a fibre mapping which is anatomically meaningful. Once again, computational anatomy appears as a key component of the clinical image analysis work-flow.

1.5 References

• [Andrade et al., 2001] Andrade A., Kherif F., Mangin J.-F, Worsley K., Paradis A.-L., Simon O., Dehaene S. and Poline J.-B.. Detection of fMRI activation using cortical surface mapping. Hum. Brain Mapp., 2001
• [Ashburner et Friston, 2000] Ashburner, J. Friston, K. J. Voxel-Based Morphometry-The Methods. NeuroImage 11(6):805-821, 2000.
• [Arsigny et al., 2006] Arsigny V., Commowick O., Pennec X., and Ayache N. A Log-Euclidean Framework for Statistics on Diffeomorphisms. In Proc. of MICCAI 2006
• [Bajcsy and Kovacic, 1989] Bajcsy R., Kovacic S., Multiresolution Elastic Matching, CVGIP, vol. 46, 1989, p . 1-21 .
• [Bondiau, 2004] Bondiau P.Y. Mise en oeuvre et évaluation d'outils de fusion d'image en radiothérapie. Thèse de sciences, Université de Nice-Sophia Antipolis, November 2004
• [Bookstein, 1978] Bookstein F.L., The Measurement of Biological Shape and Shape Change, in Lecture Notes in Biomathematics, vol. 24. New York: Springer-Verlag, 1978.
• [Collins et al, 1998] Collins D.L., Zijdenbos A.P., Kollokian V., Sled J.G., Kabani N.J., Holmes C.J., Evans A.C. Design and Construction of a Realistic Digital Brain Phantom" IEEE Transactions on Medical Imaging, vol.17, No.3, p.463--468, June 1998.
• [Christensen et al., 1997] Christensen G., Joshi S., Miller M., Volumetric transformation of brain anatomy, IEEE Trans. Med. Imag., 16, 864--877, 1997
• [Commowick 2005] Commowick O., Stefanescu R., Fillard P., Arsigny V., Ayache N., Pennec X., and Malandain G. Incorporating Statistical Measures of Anatomical Variability in Atlas-to-Subject Registration for Conformal Brain Radiotherapy. In J. Duncan and G. Gerig, editors, Proceedings of MICCAI 2005, volume 3750 of LNCS, pages 927-934, 2005. Springer Verlag
• [Commowick, 2007] Commowick O., Design and Use of Anatomical Atlases for Radiotherapy. PhD Thesis, Nice -- Sophia-Antipolis University, February 2007
• [Evans et al, 1993] Evans A. C., Collins D. L., Mills S. R., Brown E. D., Kelly R. L., Peters T. M., 3D statistical neuroanatomical models from 305 MRI volumes, Proc. IEEE-Nuclear Science Symposium and Medical Imaging Conference, 1813-1817, 1993.
• [Fillard et al, 2005] Fillard P., Arsigny V., Pennec X., Thompson P.M., and Ayache N. Extrapolation of Sparse Tensor Fields: Application to the Modeling of Brain Variability. In Gary Christensen and Milan Sonka, editors, Proc. of Information Processing in Medical Imaging 2005 (IPMI'05),2005. Springer.
• [Fillard et al, 2007] Pierre Fillard, Vincent Arsigny, Xavier Pennec, Kiralee M. Hayashi, Paul M. Thompson, and Nicholas Ayache. Measuring Brain Variability by Extrapolating Sparse Tensor Fields Measured on Sulcal Lines. Neuroimage, 34(2):639-650, January 2007. Note: Also as INRIA Research Report 5887, April 2006. PMID: 17113311.
• [Fleute and Lavallée, 1998] Fleute, M. and Lavallée S. Building a Complete Surface Model from Sparse Data Using Statistical Shape Models: Application to Computer Assisted Knee Surgery. Proc. of Medical Image Computing and Computer-Assisted Interventation (MICCAI'98), LNCS Vol. 1496, p.879-??, 1998.
• [Le Goualher et al. 1999] Le Goualher G., Procyk E., Collins L., Venegopal R., Barillot C., Evans A, Automated extraction and variability analysis of sulcal neuroanatomy. IEEE Transactions on Medical Imaging, 1999
• [Lester, 1999] Lester, H., Arridge, S.R., Jansons, K.M., Lemieux, L., Hajnal J.V. and Oatridge, A. Non-linear Registration with the Variable Viscosity Fluid Algorithm. Proc of IPMI'99, LNCS 1613, p.238-251, 1999.
• [Mangin et al., 2004] Mangin J.-F., Riviere D., Cachia A., Duchesnay E., Cointepas Y., Papadopoulos-Orfanos D., Collins D., Evans A., and Regis J. Object-based morphometry of the cerebral cortex. IEEE Transactions on Medical 23(8) 23:968-982, 2004. Imaging.
• [Mazziotta et al, 2001] Mazziotta J, Toga A, Evans A, Fox P, Lancaster J, Zilles K, Woods R, Paus T, Simpson G, Pike B, Holmes C, Collins L, Thompson P, MacDonald D, Iacoboni M, Schormann T, Amunts K, Palomero-Gallagher N, Geyer S, Parsons L, Narr K, Kabani N, Le Goualher G, Boomsma D, Cannon T, Kawashima R, Mazoyer B. 2001. A probabilistic atlas and reference system for the human brain: International Consortium for Brain Mapping (ICBM). Philos Trans R Soc Lond B Biol Sci 356:1293--1322.
• [Miller et al., 2002] Miller M., Trouvé A., Younes L., On the Metrics and Euler-Lagrange Equations of Computational Anatomy, Annual Review of biomedical Engineering, 2002
• [Modersitzki, 2004] Modersitzki J., Numerical Methods for Image Registration, Oxford University Press Series: Numerical Mathematics and Scientific Computation, 2004
• [Pennec et al, 1999] Xavier Pennec, Pascal Cachier, and Nicholas Ayache. Understanding the ``Demon's Algorithm'': 3D Non-Rigid registration by Gradient Descent. In C. Taylor and A. Colchester, editors, Proc. MICCAI'99, volume 1679 of LNCS, Cambridge, UK, pages 597-605, September 1999. Springer Verlag.
• [Pennec, 2006] Pennec X. Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements. Journal of Mathematical Imaging and Vision, 25(1):127-154, July 2006.
• [Pennec et al, 2006] Pennec X., Fillard P. and Ayache N., A Riemannian Framework for Tensor Computing. International Journal of Computer Vision, 66(1):41-66, January 2006
• [Rajamani et al, 2004] Rajamani, K.T. Joshi, S.C. Styner, M.A. Bone model morphing for enhanced surgical visualization. In Proc of IEEE Symp. on Biomedical Imaging: Nano to Macro (ISBI) 2004, Vol. 2, p.1255-1258, April 2004.
• [Rueckert et al, 2003] Rueckert D., Frangi A.F., Schnabel J.A., Automatic construction of 3-D statistical deformation models of the brain using nonrigid registration, IEEE Transactions on Medical Imaging 22,p. 1014- 1025, 2003
• [Stefanescu, 2004] Radu Stefanescu, Xavier Pennec, and Nicholas Ayache. Grid Powered Nonlinear Image Registration with Locally Adaptive Regularization. Medical Image Analysis, 8(3):325-342, September 2004. PMID: 15450226.
• [Stefanescu et al, 2005] Stefanescu R., Pennec X., and Ayache N. A Grid Service for the Interactive Use of a Parallel Non-Rigid Registration Algorithm of Medical Images. Methods of Information in Medicine, 44(2), 2005.
• [Talairach and Tournoux, 1988] Talairach J., Tournoux P. Co-Planar Stereotaxic Atlas of the Human Brain: 3-dimensional Proportional System : an Approach to Cerebral Imaging". Thieme Medical Publishers, New York, 1988.
• [Thirion, 1998] Thirion J.P. Image matching as a diffusion process: An analogy with maxwell's demon. Medical Imag. Analysis, 2:243--260, 1998.
• [Thompson et al., 1997] Thompson P.M., MacDonald D., Mega M.S., Holmes C.J., Evans A.C., Toga A.W., Detection and Mapping of Abnormal Brain Structure with a Probabilistic Atlas of Cortical Surfaces, Journal of Computer Assisted Tomography 1997
• [Thompson et al., 2000] Thompson, P., Mega, M., Narr, K., Sowell, E., Blanton, R., and Toga, A. (2000). Brain image analysis and atlas construction. Handbook of Medical Image Proc. and Analysis, chapter 17. SPIE.
• [Trouvé 1998] Trouvé A., Diffeomorphisms Groups and Pattern Matching in Image Analysis, Int. J. Computer Vision, 1998
• [Vaillant et al., 2007] Vaillant M., Qiu A., Glaunès J., Miller M., Diffeomorphic metric surface mapping in subregion of the superior temporal gyrus, NeuroImage, 2007

Computational Anatomy of the Brain

2. Morphometry of the Cortex Inferred from Sulcal Lines

2.1 A point-wise Line Correspondence Approach

In the context of the INRIA-associated team BrainAtlas with Paul Thompson at the Laboratory of Neuro Imaging, University ofCalifornia at Los Angeles (LONI), we proposed with Pierre Fillard to infer the variability of the brain from a dataset of precisely delineated anatomicalstructures (sulcal lines) on the cerebral cortex. We model the first- and second-order momentsof sulci by an average sulcal curve and a sparse field of covariance tensors along these curves. The covariance matrices are then extrapolated to the whole brainusing a harmonic diffusion PDE on the manifold of tensor fields. As a result, we obtain a dense 3D variability map, which proves to be in accordance with previouslypublished results on smaller samples of subjects. Statistical tests demonstrated that our model was globally able to recover some missing information. Preliminaryresults on the correlation between local and distant displacements indicate that the displacement of the symmetric point is correlated. Other long-distancecorrelations also appear, but their statistical significance still needs to be established.

Figure 3. From sulcal lines in a population to the brain variability; (left) sulcal lines of eighty subjects in green with the mean sulcal lines in red; (middle) variability measured along the mean sulcal lines (covariance matrix at one sigma); (right) the colour encodes the amount of variability everywhere on the cortex after the extrapolation of the variability tensors onto the whole 3D space, going from blue for the least variable places to red for the most varying ones. Images realised by P. Fillard.

2.2 A new Correspondence-less Diffeomorphic Method

Figure 4. Labelled lines registration using a model based on currents. Dark blue line is transported to the green line that matches the red line. The distance between red and green lines (i.e. the precision of the matching) is computed although the sampling of each line is different, in particular no point to point correspondences are imposed.

Secondly, the large deformation framework is rooted in the paradigm of Grenander's group action approach for modelling objects [Grenander 1993, Miller et al, 2006]. The deformations are searched in a subgroup of the group of diffeomorphisms (smooth and invertible mappings). Such deformations are defined in the whole space and act on currents. In the case of curves and surfaces this action coincides with their natural geometrical transportation.

Finally, the registration is found by minimising a cost function that makes a compromise between the regularity of the diffeomorphism and the fidelity to data computed via the distance between currents. In this setting, we do not only find an optimal matching between two lines sets but also recover a global consistent deformation of the underlying biological material on which the lines are drawn. As an example the following video show the registration of a mean brain into a subject's one constrained by the sulcal lines.

This framework has the advantage of being generative. This means that we can compute samples of sulcal lines according to the Gaussian model we have learnt on the data. In particular, we can compute a tangent-PCA on the displacement field and retrieve the different modes of deformations of the database. This gives a description of the major trends of geometrical variations within the studied population. The following video shows the first deformation mode of the database: according to the Gaussian prior, this deformation is the principal trend of variations within the population.

2.3 Comparison of Both Methods

In the correspondence-less diffeomorphic approach (called here CLD) we can take advantage of the tangent space representation of the diffeomorphisms (as in [Vaillant et al, 2004]) to compute similar maps as those obtained by the previous method set up by P. Fillard and based on point-wise lines correspondences (PwLC). These maps of variability enable to emphasise the differences between both approach. These differences are explained by the hypothesis on which each method is based.

Since the CLD's approach 'integrate' the information in a spatially consistent and smooth way, the retrieved variability is naturally more regular. This is particularly visible at lines crossing or between points that belong to different lines as highlighted in figure 5. The modelling based on the currents enables also to give an alternative to the systematically under-estimated tangential variability in the PwLC's approach as explained by the figure 6. Beyond these two obvious differences, a deeper analysis can also be performed to analyse more subtle variations as presented in [Durrleman et al, 2007].

Figure 5. At every sampling points, ellipsoids represent the covariance square root of the empirical covariance matrix. In PwLC method, the extremal points are supposed to be matched leading to a high variability at lines extremities (area 1, right). This is avoided by the approach based on currents (area 1, left). In PwLC method each line is registered individually. At line crossing, the variance can vary dramatically between two close points that belong to two different lines (area 2, right). The global regularity constraint on the left hand side leads to smoother results.

Figure 6. An alternative to the aperture problem. On the variability maps, a tangential variability is measured on the top of the brain (area 3) by the CLD method while it is intentionally minimised with PwLC. On the covariance matrices of this region, we can see that the variability is mainly longitudinal. Large tensors at lines extremities with PwLC are removed before extrapolation. The remaining tensors do not capture tangential variability and lead to small variability in this area.

2.4 References

• [Durrleman et al, 2007] Durrleman S., Pennec X., Trouvé A., Ayache N., Measuring Brain Variability via Sulcal Lines Registration: A Diffeomorphic Approach, MICCAI 2007, to appear.
• [Grenander, 1993] Grenander U., General Pattern Theory, Oxford Science Publications, 1993
• [Miller et al, 2006] Miller M., Trouvé A., Younes L., Geodesic Shooting for Computational Anatomy. Journal of Mathematical Imaging and Vision, 2006.
• [Pennec, et al 2005] Xavier Pennec, Radu Stefanescu, Vincent Arsigny, Pierre Fillard, and Nicholas Ayache. Riemannian Elasticity: A statistical regularization framework for non-linear registration. In J. Duncan and G. Gerig, editors, Proceedings of the 8th Int. Conf. on Medical Image Computing and Computer-Assisted Intervention - MICCAI 2005, Part II, volume 3750 of LNCS, Palm Springs, CA, USA, October 26-29, pages 943-950, 2005. Springer Verlag. PMID: 16686051.
• [Pennec, 2006] Xavier Pennec. Left-Invariant Riemannian Elasticity: a distance on Shape Diffeomorphisms?. In X. Pennec and S. Joshi, editors, Proc. of the International Workshop on the Mathematical Foundations of Computational Anatomy (MFCA-2006), pages 1-13, 2006.
• [Vaillant et al, 2004] Vaillant M., Miller M., Younes L., Trouv�A., Statistics on diffeomorphisms via tangent space representation, NeuroImage 23, 161-169, 2004
• [Vaillant and Glaunès, 2005] Vaillant M., Glaunès J., Surface Matching via Currents, Proceedings of Information Processing in Medical Imaging. Volume 3565 of LNCS. Springer (2005)

Computational Anatomy of the Brain

3. A Statistical Framework for Anatomical Curves and Surfaces

3.1 Specific Approximation Tools for Groupwise Statistics

As emphasized in Section 1, one the major goal of Computational Anatomy is to compute, process and visualize statistics on large database of anatomical curves or surfaces. As shown on the sulcal lines in Section 2, modeling such geometrical primitives with currents avoids to rely on point-correspondence method between structures, which may introduce bias in the statistical estimation. It also allows to use all the information of the curve or surface without having to chose a priori a reduced set of characteristic measurements (length, area, volume, curvature, etc.).The framework based on currents has been proved to be powerful for pairwise registration of brain surfaces. However, new numerical schemes are required to process groupwise statistics due to an increasing complexity when the size of the database is growing. For example, the registration scheme for currents has a polynomial complexity in the number of points in shapes. In the space of current the mean or the first mode contains as many points as the total number of points within the dataset. Registering the mean shape therefore becomes untractable as the number of instances in the database becomes important (which is what we aim at for good statistics).

To solve that problem, we propose an efficient computational framework that overcomes this limitation by providing a sparse representation of currents at any desired accuracy [Durrleman et al., 2008b]. Actually, the statistics of currents such as mean or principal modes often have a highly redundant representation. Our algorithm builds on ideas from the approximation theory previously developed to decompose images in wavelet bases [Mallat and Zhang, 1993]. To the very best of our knowledge, this is the first time that these signal processing techniques are used in geometric shape analysis. As a result, the method finds an adapted basis, which minimizes in some sense the redundancy of the representation. Besides the computational improvement, this sparse representation offers a way to visualize and interpret statistics on currents. Following experimental results clearly demonstrate the interest of our method.

3.2 Groupwise Statistics on Sets of Curves

The sulci are the fissures on the brain surface and they are often used to measure anatomical differences between subjects, as emphasized in Section 2. We use here a set of 70 sulci delineated in 34 subjects. For each sulcal line, we give a sparse representation of the mean current with an approximation error below 5% of the variance. Results are shown in Fig. 1 for all 70 sulci. A close-up of the Sylvian Fissure of the right hemisphere is shown in Fig. 2. The initial mean fissure contains 899 segments (i.e. the number of segments of all lines), whereas the final approximation requires only 54 momenta. This gives a compression ratio of 94%. Considering all sulci, the compression ratio is on average: 94.8%. We see also that our mean is in good agreement with other mean curves computing from B-spline representation [Fillard et al., 2007].

Then, we compute a PCA on the lines set to retrieve its first eigenmodes. The eigenmodes are given as linear combinations of the input currents. We approximate these combinations with an error below 5% of the variance. The first eigenmodes of the Sylvian Fissure of the right hemisphere is shown in Fig. 2. This mode captures mainly the spreading of the lines set.

Figure 1. Mean Curves from Sulcal Lines Original data: 70 sulcal lines delineated in 34 subject (left). Mean curves computed by our method in red and by a point-wise line correspondence method in blue (right). Both means are in good agreement.

Figure 2. Close-up on the Sylvian Fissure (right hemisphere) The original 34 curves (left), our mean current in red (middle) and the first mode of variations in green (right). The latter shows mainly the spreading of the lines set. The mean line computed from a point-wise line correspondence approach is superimposed in blue.

3.3 Groupwise Statistics on Sets of Surfaces

We compute the mean current of 10 meshes of deep brain structures: Caudate, Putamen, Globus Pallidus, Amygdala and Hippocampus for each hemisphere for 50 autistics children of age 2 [Hazlett et al., 2005]. For surfaces, we represent a momenta by a small triangle whose normal is the momenta. Results are shown in Fig. 3. The compression ratio achieved for these 10 surfaces is on average of 99.96%. Representing the mean requires 1.2 Kb in our framework, versus 8 Mb originally. Fig. 3 (right panel) shows that the quality of approximation remains good until very high compression ratio. The deformation of a mean obtained from only 3 shape instances, each with 15 000 points, which was previously taking 10 hours, is now taking about 5 minutes (using the same code as in [Vaillant and Glaunes, 2005; Durrleman et al., 2008a]. For the full set of 50 instances, deforming the former still requires 5 minutes while it is not feasible to deform the later without high performance computing.

Similarly, we compute mean from meshes of 7 controls. The difference between both means is a current that we approximate: the arrows of Fig.4 are the first 10 estimated momenta of this difference, suggesting that the autistic mean is more curved at the Hippocampus' extremity and thicker in the middle. Such results still need to be confirmed by rigorous statistical tests.

Figure 3. Mean of subcortical structures from 50 autistics children of age 2 Left: One instance of the database (left): Caudate (dark blue), Putamen (light blue), Globus Pallidus (green), Amygdala (yellow) and Hippocampus (red). Middle: the mean of 50 autistics, approximated with an approximation error below 5% of the variance. Right: the diagram shows that very high compression ratio of the mean can be achieved, while the quality of the approximation remains very good. The red cross corresponds to compression ratio of 5% of the variance, shown in the middle panel. Data courtesy of Guido Gerig and Martin Styner [Hazlett et al., 2005]

Figure 4. Difference between autistic and control mean anatomy The mean of 7 controls is computed (left). The difference between the mean of autistics and mean of control is approximated by our method, leading to the red arrows in the right part of the figure. These arrows are superimposed with the Hippocampus of a control. They show that the mean from autistics is more curved at hippocampus' extremity (area 1) and thicker in area 2.

3.4 References

• [Durrleman et al, 2008a] Durrleman S., Pennec X., Trouvé A., Thompson P.M., Ayache N., Inferring brain variability from diffeomorphic deformations of currents: an integrative approach. Medical Image Analysis, 12 (5):626-637, 2008
• [Durrleman et al, 2008b] Durrleman S., Pennec X., Trouvé A., Ayache N., Sparse Approximation of Currents for Statistics on Curves and Surfaces. In Proceedings of Medical Image Computing and Computer Assisted Intervention (MICCAI), Part II, volume 5242 of LNCS, New-York, USA, pages 390-398, Springer (2008).
• [Fillard et al., 2007] Fillard P., Arsigny V., Pennec X., Hayashi K., Thompson P.M., Ayache N., Measuring Brain Variability by extrapolating sparse tensor fields measured on sulcal lines. NeuroImage 34 (2) pages 639-650, 2007
• [Hazlett et al., 2005] Hazlett H., Poe M., Gerig G., Smith R., Provenzale J., Ross A., Gilmore J., Piven J., Magnetic Resonance Imaging and Head Circumference Study of Brain Size in Autism, The Archives of General Psychiatry 62, pages 1366-1376, 2005
• [Mallat and Zhang, 1993] Mallat S., Zhang Z., Matching Pursuits with time-frequency dictionaries, IEEE Transactions on Signal Processing 41 (12), pages 3397-3415, 1993
• [Vaillant and Glaunès, 2005] Vaillant M., Glaunès J., Surface Matching via Currents, Proceedings of Information Processing in Medical Imaging. Volume 3565 of LNCS. Springer (2005)

Computational Anatomy of the Brain

4. A Generative Model of Brain Populations Variability

4.1 A glimpse at the existing methods

To study the variability of the brain in populations, a first group of studies focused on voxel-based statistics [Hermoye et al. 2006, Giedd et al. 1999, Thompson et al. 2001]. In this case, one segments the brain into several regions of interest (grey matter, white matter, cerebro-spinal fluid, corpus callosum, thalamus) and one look at the evolution of the volume of such predefined regions of interest. In [Gogtay et al 2004], local density of grey matter is measured on each point of the cortex: this analysis at the voxel's level avoids the a priori definition of regions of interest that appears afterwards in the density maps. Some studies use DTI to perform longitudinal analysis like in [Barnea-Goraly et al., 2005] for example. However, the analysis is also mainly limited to the evolution of a single value: the fractional anisotropy calculated on anatomical landmarks. In every case, statistics on the evolution are performed on the basis of scatter plots. One computes linear or quadratic regression and then performs statistical tests (like Hotteling one) on the retrieved parameters to emphasize possible significant differences in the growth according to external features like gender, handedness, and ethnicity, or according to some developmental disorders like Childhood-onset schizophrenia or Williams syndrome. Eventually, like in [Gilmore et al, 2007], one measures anatomical differences across the population thanks to statistical tests performed on geometrical features such as cerebral asymmetry.

Another approach consists in deformation-based statistics ([Ashburner et al. 1998]) as initiated in [Chung et al. 2000a, Chung et al. 2000b, and Chung et al.2003] for temporal evolutions. In this approach the different acquisitions of the same subject are registered into the baseline leading to continuously time-dependent displacement field. Given this field for several subjects, one can fit a spatially varying linear deformation rate. Then, Hotelling's tests on the residuals determines the regions of brain where growth or atrophy is statistically significant. Notice that [Chung et al. 2003] synthesizes the evolution of surfaces by the evolution of geometrical invariants: area, volume and curvature. As far as we know, this is the first attempt to describe the anatomical variability in a more geometrical sense and to build a joint model that takes into account the variability both across a population and within a given subject's growth.

All these recent approaches enable to achieve some promising results. The specific role of the grey-matter evolution was highlighted and maturation of the brain was quantified. In particular, the process of myelination and its consequences was broadly discussed. Moreover, one emphasizes some major trends in the brain maturation and statistical tests enable to show significant differences between healthy and abnormal maturation. Such tests also reveal that external features like handedness have a significant impact on brain growth. However, the description of the brain maturation and the differences between different evolutions remain very simple. In most cases the whole geometry of the brain was summarized by the measure of its volume. Actually, the maturation issue was considered so far only as a problem of brain's volume growth. One would like to go one step further and describe the brain's anatomy more accurately. Looking at the volume evolution of sub-regions of the shape could not be satisfactory: we need continuous models of shapes and of their evolutions in order to provide a non-supervised geometric description of the difference between the evolution of two shapes. The central issue is actually to describe how the brain is growing, how anatomical structures like sulcal landmarks or cortical surfaces are developing through the time, how they differ one from another in a geometrical point of view.

A second feature that was not yet taken into account is the joint cross-sectional and longitudinal variability. Age is mainly considered as an absolute variable and data of different subjects are matched according to the age. The question of a possibly developmental delay or advance was not put so far.

Eventually, most of the models proposed in the literature are descriptive. The approaches based on statistical tests only state that there is a difference at the population level without providing a statistical model that explains the observed data. We believe that the conclusions would be more convincing if one could generate new observations according to a learnt statistical model. This would make obvious the variability captured by the model. Moreover such generative models could be also predictive: the classification parameters determined on a population in the learning phase can be turned into a classifiers on individual cases that could be used for diagnosis.

Our goal is therefore to define theoretical tools addressing the aforementioned limitations of the current work and to define statistical models that could be calibrated and tested on existing databases. Our strategy is to address the problem in two steps. First, we will define a generative statistical model for subjects of the same age. This model will able to capture a large range of geometrical variations, without reducing the shapes as a sets of arbitrary features. It will also combine deformation-based statistics and voxel-based statistics. Once this model will be defined and evaluated on real data, it will be extended in the next section to account for temporal evolutions. The resulting framework will therefore address both the geometrical and the temporal aspects.

4.2 A Generative Statistical Model at a Given Time-Point

In the medical imaging field, frameworks have been proposed to build atlases from databases of medical images [Joshi et al., 2004; Avants and Gee, 2004; Zollei et al., 2005] or from anatomical curves or surfaces [Glaunes and Joshi, 2006; Chui et al., 2004]. In both cases, the underlying idea is to estimate a "mean anatomy" (called also template or atlas in the literature), and then learn how this mean model deforms within a given population. We choose here to base our statistical estimation on a forward model, as pioneered in [Allassonniere et al., 2007; Ma et al., 2008]. This model appears to be particularly well suited for comparison between new data and the atlas, as explained in [Durrleman et al., 2008]. In this setting the observations are considered as a random perturbation of an unknown template, plus a residual perturbation.

The observations are given as discrete sampled objects. The template models an average continuous "ideal" biological material. The deformations of the template within the population capture the geometrical variability such as torque, stretching or shrinking effects. The residual perturbations is what remain once the geometrical variability has been discounted. It is measured by the differences between the observations and the deformed template (i.e. what cannot be captured by a regular deformation). It captures the variability in terms of "texture", which contains, for instance, change of topology, material creation or deletion, physical or numerical noise.

To estimate consistently the template, the deformations of the template in the population and the residual perturbation, we propose a Maximum A Posteriori approach. This approach takes advantage of the correspondence-less diffeomorphic approach detailed in Section 2 to deform the template to each subject's anatomy. Figure 1 present a template of white matter fiber bundles segmented in 6 subjects using MedINRIA.

Figure 1: Template of five fibers bundles segmented in 6 subjects Left: One subject. the corticospinal tract (blue), the corticobulbar tract (yellow), the callosal fibers (red), the left and right arcuate fasciculi (green). The frontal (middle) and lateral (right) views of the template learned from the dataset. This estimation relies on our modeling based on currents: the comparison between this template and each data do not assume either point or fiber correspondences between subjects. It is robust to fiber interruption. (data courtesy of D. Ducreux, MD, Bicetre Hospital, Paris. This work is a collaboration with P.Fillard, CEA, Neurospin.)

Once, the atlas is build, we can describe the variability within the population by analyzing (1) the geometrical variability captured by the deformations and (2) the variability in terms of texture captured by the residual perturbation. This is done by performing PCA on the deformations parameters and the residual currents respectively. This last PCA in the space of shapes is performed as shown in Section 3. Such a statistical analysis has been performed on our dataset of fiber bundles and results are shown in Fig. 2 and Fig. 3.

Figure 2: First deformation mode of the corticobulbar bundle (view from top). This mode captures the diffeomorphic variability of the population around the prototype bundle (middle) at plus sigma (right) and minus sigma (left). It shows a torque at the basis of the bundle and a stretching/shrinking effect of the left and right parts of the bundle. Torque may be related to the well-known brain torque.

Figure 3: First "texture" mode of the corticobulbar bundle. This mode captures the residual variability of the population, once the diffeomorphic variability has been discounted. This mode is added to (right) or removed from (left) the prototype bundle (middle). The mode at plus (resp. minus) sigma shows that the left (resp. right) part of the bundle becomes thicker, while its right (resp. left) part becomes thiner. Further investigation should determine whether this asymmetry is a true anatomical feature or an artifact of the tracking method.

We can also use the atlas to generate new data according to the law we learn on data. We simulate an instance of the deformation, we deform the template according to this deformation and we add an instance of the residual to the deformed template. This is illustrated in Fig. 4.

Figure 4: Simulation of sythetic corticobulbar bundle Our statistical model learn how to simulate new bundles from the prototype. First we deform the template (left) according to the first mode of deformation (middle). We then add the first texture mode to this deformed template, to give a new instance of the bundle (right). We can compare this simulation with the real data to better understand how our model interpolates the variability in-between the train data. This illustrates the generative property of our modeling.

We can compute how likely a new observation could be according to our atlas. Given a new observation, we register the template to this data. This decomposes the variability between the observation and the template into a geometrical deformation and a residual perturbation. Then, we can measure the likelihood of the deformation's parameter and the residual current according to the Gaussian law we learnt on the parameter space and the space of currents. This could be used to classify new subjects according to their gender, handedness, pathologies, etc., as well as to detect abnormalities via deviations from the normal variability of the organs. One can also think of using the atlas to constraint the segmentation of the organs in images of new subjects.

4.3 References

• [Allassoniere et al., 2007] Allassonniere S., Amit Y., Trouvé A., Toward a coherent statistical framework for dense deformable template estimation. Journal of the Royal Statistical Society Series B 69 (1), p. 3-29, 2007
• [Ashburner et al., 1998] Ashburner J., Hutton C., Frackowiak R., Johnsrude I., Price C., Friston K., Identifying global anatomical differences: deformation-based morphometry, Human Brain Mapping 6, 1998
• [Avants and Gee, 2004] Avants B. and Gee J., Geodesic estimation for large deformation anatomical shape averaging and interpolation, Neuroimage 23, pp. 139-150, 2004
• [Barnea-Goraly et al., 2005] Barnea-Goraly N., Menon V., Eckert M., Tamm L., Bammer R., Karchemskiy A., Dant C.C., Reiss A.L., White Matter Development During Childhood and Adolescence: A Cross-sectional Diffusion Tensor Imaging Study, Cerebral Cortex, 2005
• [Carissa et al., 2007] Carissa J.C., Gerig G., Piven J., Diffusion Tensor Imaging: Application to the Study of the Developing Brain, Journal of the American Academy of Child and AdolescentPsychiatry, 2007
• [Chui et al., 2004] Chui H., Rangarajan A., Zhang J., Leonard C., Unsupervised learning of an atlas from unlabeled point-sets, IEEE Transactions on Pattern Analysis and Machine Intelligence 26 (2), pp. 160-172, 2004
• [Chung et al., 2000a] Chung M.K., Worsley K.J., Paus T., Cherif C., Collins D.L., Giedd J.N., Rapoport J.L. and Evans A.L., A Unified Approach to Deformation-Based Morphometry, NeuroImage, 2000
• [Chung et al., 2000b] Chung M.K., Worsley K.J., Paus T., Cherif C., Collins D.L., Rapoport J.L. and Evans A.L., Statistical Analysis of Local Volume Change, with an Application to Brain Growth, NeuroImage 2000
• [Chung et al., 2003] Chung M.K., Worsley K.J., Robbins S., Paus T., Taylor J., Giedd J.N., Rapoport J.L. and Evans A.L., Deformation-Based Surface Morphometry Applied to Gray Matter Deformation, NeuroImage:18(2), 2003
• [Durrleman et al., 2007] Durrleman S., Pennec X., Trouvé A., Ayache N., Measuring Brain Variability via Sulcal Lines Registration: A Diffeomorphic Approach, Proc. Medical Image Computing and Computer Assisted Intervention (MICCAI), LNCS 4791, pp. 675-682, 2007
• [Durrleman et al., 2008] Durrleman S., Pennec X., Trouvé A., Ayache N., A forward model to build unbiased atlases from curves and surfaces, Proc. of Mathematical fundations of Computational Anatomy, 2008
• [Evans et al, 2006] Evans A.L. and the Brain Development Cooperative Group, The NIH MRI study of normal brain development, NeuroImage 30(1):184-202, March 2006
• [Gerig et al, 2006] Gerig G., Joshi S., Fletcher T., Gorczowski K., Xu S., Pizer S.M. and Styner M. Statistics of population of images and its embedded objects: driving applications in neuroimaging. In Proc. of IEEE ISBI, pp1120-1123, April 2006.
• [Giedd et al., 1999] Giedd J.N., Blumenthal J., Jeffries N.O., Rajapakse J.C., Vaituzis C., Liu H., Berry Y.C., Tobin M., Nelson J., Castellanos X., Development of the human corpus callosum during chilhood and adolescence: a longitudinal MRI study, Prog. Neuro-Psychopharmacol. & Biol. Psychit. 23, 1999
• [Gilmore et al., 2007] Gilmore J., Lin W., Prastawa M., Looney C., Vetsa Y.S., Evans D., Smith J., Hamer R., Lieberman J., Knickmeyer R., Guido G., Regional Gray Matter Growth, Sexual Dimorphism and Cerebral Asymmetry in the Neonatal Brain, The Journal of Neuroscience 27, 2007
• [Glaunès and Joshi, 2006] Glaunès J. and Joshi S., Template estimation from unlabeled point set data and surfaces for Computational Anatomy, Proc. of Mathematical Foundations of Computational Anatomy, 2006
• [Gogtay et al., 2004] Gogtay N., Giedd J.N., Lusk L., Hayashi K., Greenstein D., Vaituzis C., Nugent III T.F., Herman D.H., Clasen L.S, Toga A., Rapoport J.L., Thompson P.M., Dynamic mapping of human cortical development during childhood through early childhood, P.N.A.S., 2004
• [Hazlett et al, 2005] Hazlett, H.C., Poe, M., gerig, G., Gimpel Smith R.,Provenzale J., Ross Å, Gilmore J.H. and Piven J. Magnetic resonnance imaging and head circumference study of brain size in autism. In The Archives of General Psychiatry, vol. 62, pp.1366-1376, December 2005.
• [Hermoye et al., 2006] Hermoye L., Saint-Martin C., Cosnard G., Lee, S.-K., Kim J., Nasogne M.-C., Menten R., Clapuyt P., Donohue P.K., Hua K., Wakana S., Jiang H., van Zijl P.C.M., Mori S., Pediatric diffusion tensor imaging: Normal database and observation of the white matter maturation in early childhood, NeuroImage 29, 2006
• [Joshi et al., 2004] Joshi S., Davis B., Jomier M., Gerig G., Unbiased diffeomorphic atlas construction for Computational Anatomy, Neuroimage 23, pp. 151-160, 2004
• [Ma et al, 2008] Ma J., Miller M., Trouvé A., Younes L., Bayesian template estimation in Computational Anatomy, Neuroimage 42 (1), pp. 252-261, 2008
• [Paus et al., 2001] Paus T., Collins D.L., Evans A.C., Leonard G., Pike B., Zijdenbos A., Maturation of white matter in the human brain: A review of magnetic resonance studies, Brain Research Bulletin, 54, 2001
• [Prastawa et al., 2005] Prastawa M., Gilmore J.H., Lin W., Gerig G., Automatic segmentation of MR images of the developing newborn brain, Medical Image Analysis: 9, 2005
• [Thompson et al., 2001] Thompson P.M., Vidal C., Giedd J.N., Gochman P., Blumenthal J., Nicolson R., Toga A.W., Rappoport J.L., Mapping adolescent brain change reveals dynamic wave of accelerated gray matter loss in very early-onset schizophrenia, P.N.A.S. 98, 2001
• [Vaillant and Glaunes, 2005] Vaillant M., Glaunes J., Surface Matching via Currents, Proc. of Information Processing in Medical Imaging, 2005
• [Zollei et al., 2005] Zollei L., Learned-Miller E., Grimson E., Wells W., Efficient population registration of data, Proc. of Computer Vision for Biomedical Image Applications, LNCS 3765, pp. 291-301, 2005

Computational Anatomy of the Brain

5. Spatiotemporal Atlas Estimation in Longitudinal Datasets

In the context of the Health-e-Child project, we aim at understanding the child brain maturation from an anatomical point of view. As emphasized by [Toga et al., 2006], "the dynamic course of brain maturation is one of the most fascinating aspects of the human condition. Although brain change and adaptation are part of a lifelong process, the earliest phases of maturation -during foetal development and childhood- are perhaps the most dramatic and important. Indeed, much of the potential and many of the vulnerabilities of the brain might, in part, depend on the first two decades of life." The newborn brain will grow up to a factor 3 or 4. During this period, the axons will be progressively myelinated and the brain will specialize according to both its genetic program and its environmental influences: the number of synaptic connections will dramatically increase whereas some connections will be removed to give a 'more efficient' set of connections.

Whereas the detailed maturation of the brain is obviously specific to every child we believe that there are some common features and that the normal maturation follows the same underlying standard process. Our goal is therefore to describe such a mean scenario and to analyze the deviations of a given brain maturation from this mean. With such statistical models we would be able to divide a given population into statistically consistent clusters, to classify new observations according to models learnt on large datasets, to analyze the correlations between brain growth and extrinsic features like gender or handedness, to detect abnormal evolutions as deviations from the normal models, to register atlas to patient's specific anatomy. Here, we face not only a high shape variability among different subjects (inter-subject variability) but also high variations of the same brain trough the time. The central issue is therefore to model the variability of geometric features evolutions (in a space of (3+1) dimensions) instead of the variability of the geometric feature at a given age (in a 3D space). Indeed, it seems not to be relevant to integrate individual measures of child brains independently at each age since the brain maturation stages appear at different ages, with different speed for every child. Instead of computing a mean shape of a given anatomical structure and analyzing the deviation of the anatomy around this mean, we would better compute a mean scenario and analyze the deviation both in space and time of the maturation process around this mean scenario.

In the previous section, we develop a complete statistical model that decomposes the variability of anatomical structure into a geometrical and a texture part. All the available information is taken into account: we can simulate new data according to our model. The experiments on the fiber bundles shows undoubtedly that both parts contains interesting anatomical features. They show the large range of possible variations that our model is able to retrieve. It comes from the fact that our model is not based on strong assumptions on the nature of the variability we are looking for.

However, this analysis could not be performed directly on temporal data that have not been registered in time. Our current work focus on extending this framework to account for temporal data. An observation is a couple shape plus time. Our extended framework will consist of a mean scenario (an evolving template) which account for all the temporal evolution and spatial deformation of this template which accounts for the inter-subject variability at every time points. Besides the usual spatial variability across subjects, the framework also enables to describe the variability in terms of speed growth via temporal registration. The whole section is based on the work presented in [Durrleman et al. MICCAI 2009].

5.1 Methodology for Statistics on Longitudinal Data

Many frameworks has been already proposed in medical imaging to analyze the anatomical variability of 3D structures like images, curves or surfaces. Less attention has been paid to the variability of longitudinal data (several subjects scanned several times). In [Qiu et al. 2008], the evolution between two shapes is modeled by a geodesic deformation, which cannot be used for more than two data per subjects. In [Gorczowski et al 2009], shape growth is measured via the evolution of extracted features like volumes, shape or pose parameters. In [Davis et al, 2007] and [Khan & Beg 2008], a temporal regression is proposed globally for a population, but this does not allow inter-subject comparisons. In cardiac motion analysis [Chandrashekara et al. 2003, Peyrat et al. 2008], spatiotemporal registration relies on 3D-registrations between images of the same moment of the cardiac cycle and between two consecutive time-points. These works rely on time-point correspondence and do not call the labels of the time-points into question. By contrast, in longitudinal studies, subjects are scanned at ages which do not necessarily correspond. Moreover, evolutions may be delayed or advanced within a population, a key feature that we precisely aim at detecting. In [Declerck et al. 1998, Perperidis et al. 2005], deformation of cardiac motion are proposed both in space and time but they require a fine temporal sampling of the motion, whereas only few acquisitions per subjects are available in most longitudinal studies. In this work, we propose to use a regression model to estimate a continuous evolution from data sparsely distributed in time and spatiotemporal deformations which register jointly both the 3D geometry and the scenario of evolution. Geometrical data are modeled as currents to avoid assuming point correspondence between structures. Large deformations are used which gives a rigorous framework for statistics on deformations and atlas construction. From longitudinal data, we estimate consistently the most likely scenario of evolution and its spatiotemporal variability within the population.

Here, we call longitudinal data a set of geometrical data (curves or surfaces, called here shapes), acquired from different subjects scanned at several time-points. We assume that the successive data of a given subject are temporal samples of a continuous evolution. We propose therefore a regression model which computes a continuous evolution which matches the data of the subject at the corresponding time points (Fig. 1). This continuous evolution allows us to compare two subjects at a given age, even if one subject has not been scanned at this age.We can also analyze how the shape varies near this age to detect possible developmental delays. We define then the spatiotemporal deformation of a continuous evolution, which consists of two deformations: (1) a morphological deformation (of the 3D space) which changes the geometry of every frame of the evolution independently of the time point and (2) a time change function (deformation of the time interval) which changes the dynamics of the evolution without changing the geometry of shapes. To avoid time-reversal, the time change function must be smooth and order preserving: it is a diffeomorphism of the time interval of interest. A 4D registration between two subjects looks for the most regular spatiotemporal deformation, such that the deformation of the continuous evolution inferred from the first subject maps the successive target data (Fig. 2). Eventually, we use this 4D registration framework to estimate a spatiotemporal atlas from a population, based on an 4D extension of the statistical model of [Durrleman et al. MEDIA 2009]. We look for a template and a continuous evolution of this template (called mean scenario of evolution), so that data of each subject are temporal samples of a spatiotemporal deformation of the mean scenario. A Maximum A Posteriori estimation enables to estimate consistently the template, the mean scenario and the spatiotemporal deformations of this mean scenario to each subject.

5.2 Regression Model for Shape Evolution

We want to fit a continuous shape evolution to a set of shapes (S_i) of the same subject acquired at different time points (t_i). Without loss of generality, we can assume that t_min = 0 and t_max = T. This evolving shape is equal to the baseline M₀ at time t = 0, which may be the earliest shape of this subject or a template as in Sec. 5.4. The evolution has the form: M_t = χ_t(M₀) where t varies continuously in the time interval [0, T]. For each t, χ_t is a diffeomorphism of the 3D space, such that χ₀ = Id (which leads to χ₀(M₀) = M₀). The regression (M_t) must match the observation S_i at the time-points t_i, while a rigidity constraint controls the regularity of the regression. This is achieved by minimizing:

5.3 Spatiotemporal Pairwise Registration

We assume now that we have successive shapes for the source subject S_i=S(t_i) and for the target T_j=T(t_j). As in Sec 5.2, we perform a regression on the source shapes which leads to a continuous evolution S(t). Our goal is to find a diffeomorphism of the 3D space φ and a diffeomorphism of the time-interval ψ which deform the source evolution S(t) into S'(t) = φ(S(ψ(t))) such that S'(t_j) match T(t_j). Thanks to the regression function, no correspondence is needed between the time points t_i and t_j . Formally, we minimize:

The spatial (resp. temporal) deformation φ (resp. ψ) is solution at parameter u = 1 of the flow equation ∂_u φ_u(t)= v_u^ψ(ψ(t)). The norm of the speed vector fields v^φ and v^ψ integrated between 0 and 1 defines the regularity terms Reg(φ) and Reg(ψ) respectively. As in Sec 5.2, the geometrical (resp. temporal) deformation is parametrized by momenta α (resp. β) at the points of S(o) (resp. at the poinst of T(0)), which are used as variables for the gradient descent. Details are given in [Durrleman et al. MICCAI 2009]. We used here the fact that the shapes are embedded with a vector space (the space of currents) provided with an inner product. Note that we minimize J with respect to the geometrical and the temporal parameters jointly. We do not performed alternated minimization.

5.4 Spatiotemporal Atlas Construction

We assume now that we have a set of N subjects (Sⁱ), provided each with temporal observations (Sⁱ(tⁱ_j). We are looking for a template M₀ and a mean scenario of evolution of this template M(t) = χ_t(M₀), such that the observations correspond to particular moments of a spatiotemporal deformation of the mean scenario. A Maximum A Posteriori estimation in the same setting as in [Durrleman et al. MEDIA 2009], leads to the minimization of

5.5 Numerical Experiments

Evolution of 2D Curves

This experiments relates to Sec. 5.2 and 5.3. We have five 2D-profiles of hominids skulls which consist of six lines each (source: www.bordalierinstitute.com). Our regression framework infers a continuous evolution from the Australopithecus to the Homo sapiens sapiens which matches the intermediate stages of evolution in Fig. 1.

Get the Flash Player to see this player. Video: Evolution of the skull profile of hominids. Copyright INRIA 2009.

Figure 1: Skull profile of five hominids (in red). The regression model estimates a continuous evolution (in blue) of the Australopithecus, which closely matches the data.

Then, we register the evolution {Homo habilis-Homo erectus-Homo neandertalensis} to the evolution {Homo erectus-Homo sapiens sapiens} in Fig 2.

Figure 2: Registration of the evolution {Homo habilis-erectus-neandertalensis} (in red) to the evolution {Homo erectus-sapiens sapiens} (in green), shifted to start at the same time. Top row: Regression of the source data (red) gives the continuous evolution in blue. Middle row: The geometrical part φ is applied to each blue frame. This shows morphological changes: the skull is larger, rounder and the jaw less prominent. Bottom row: The time change function ψ is applied to the evolution of the second row. The blue shapes are moved along the time axis (as shown by dashed black lines), but they are not deformed. Black arrows show that a better alignment is achieved when one accounts both for morphological changes and a change of the evolution speed.

The geometrical deformation shows that during the later evolution the jaw was less prominent and the skull larger and rounder than during the earlier evolution. The time change function shows that the later evolution occurs at a speed 1.66 times faster than the earlier evolution. This value is compatible with the growth speed of the skull during these periods (See Fig. 3): between Homo erectus and sapiens the skull volume growths at (1500 − 900)/0.7 = 860cm3/106years, whereas between Homo habilis and neandertalensis, it growths at (1500 − 600)/1.7 = 530cm3/106years, namely 1.62 times faster.

Figure 3: a- time change function ψ(t) of the registration in Fig. 2 (in black the reference ψ(t) = t). The slope of the curve measures an acceleration between evolutions, which is compatible with the growth of skull volume in b (source: www.bordalierinstitute.com).

Evolution of 3D surfaces

This experiment illustrates the complete spatiotemporal atlas construction detailed in Sec. 5.4 with surfaces extracted from brain images. We use here meshes of amygdala of the right hemisphere from 4 autistics, 4 developmental delay and 4 control children scanned twice [Hazlett et al. 2005]. Age distribution is shown in Fig 5-a. From these data registered rigidly, we infer a template, a mean scenario of evolution of this template and the spatiotemporal evolution of this mean scenario to each subject. In the setting of [Miller et al. 2002, Durrleman et al. MEDIA 2008], the diffeomorphisms are controlled by the standard deviation of Gaussian kernel set to 15 mm for χ, φand 1 year for ψ. The typical scale on currents is set to 3 mm. The trade-off parameters were set to γ_χ= 10^-3 and γ_φ=10^-6. An amygdala is typically 10 mm large. The discrete time step is set to 0.2 years.

Get the Flash Player to see this player. Video: Mean Scenario of the right Amygdala growth. Copyright INRIA 2009.

Figure 4: Mean Scenario of the right Amygdala (right lateral part). Arrows measures the differences between age t+0.2 and age t in the space of currents as in [Durrleman et al. MICCAI 2008]. From age 2 to 2.8, the evolution is mainly a torque at the posterior part; then the structure becomes thicker, mostly at the superior part between age 2.8 and 4 and at the inferior between age 4 and 6; from age 6 the evolution is a mainly a torque at the anterior part.

By inspection of the companion movie of the mean scenario, one distinguishes 4 phases during growth (See also Fig. 4). Preliminary tests do not show correlations between the morphological deformations and the pathology. From the time change functions shown in Fig. 5, we cannot conclude that a subject with pathology is systematically delayed or advanced compared to controls, even at a given age. However, the curves show that the growth speed seems to follow the same pattern, mainly an acceleration between age 2.5 and 3.5 for the autistics and between age 4 and age 5 for controls. The developmental delay also have such pattern but it occurs at a very variable age. These results suggest that the discriminative information between classes might not be inferred from the anatomical variability at a given age, but rather from variations of the growth process. These results, however, must be strengthen using larger database. The more time-points per subjects, the more constrained the mean scenario estimation. The more subjects, the more robust the statistics.

Figure 5: Temporal deformation of the mean scenario. Left: distribution of original (top) and registered (bottom) ages. Middle: time change functions for the 12 subjects. Right: First mode of variation at ± of the time change functions for each class. Autistics and controls show the same evolution pattern, but shifted in time.

5.6 Discussion and Conclusion

In this work, we proposed a generic framework to analyze variability of longitudinal data. A regression model fits a continuous evolution to successive data of one subject. 4D registrations decompose the difference between two sets of longitudinal data into a geometrical deformation and a change of the dynamics of evolution. The more acquisitions per subjects, the more constrained this decomposition. However, no constraint is imposed in terms of number and correspondence of measurement points across subjects. These pairwise registrations are used for group-wise statistics: ones estimates consistently a template, the mean evolution of this template and the spatiotemporal variability of this evolution in the population. Then, statistical measures can be derived, like the first mode of temporal deformation in Fig. 5. These first results suggest that pathologies might be characterized more by a particular scenario of evolution than by the anatomy at a given age. Our methodology can be used therefore to drive the search of new anatomical knowledge and to give characterization of pathologies in terms of organ growth scenario. This may be applied to the study of degenerative diseases or cardiac motion disorders.

5.7 References

• [Chandrashekara et al. 2003] Chandrashekara, R., Rao, A., Sanchez-Ortiz, G., Mohiaddin, R.H., Rueckert, D.: Construction of a statistical model for cardiac motion analysis using nonrigid image registration. In: Proc. of IPMI. Volume 2732 of LNCS., Springer (2003) 599-610
• [Davis et al. 2007] Davis, B., Fletcher, P., Bullitt, E., Joshi, S.: Population shape regression from random design data. In: Proc. of ICCV'07. 1-7
• [Declerck et al. 1998] Declerck, J., Feldman, J., Ayache, N.: Definition of a 4D continuous planispheric transformation for the tracking and the analysis of LV motion. Medical Image Analysis 4(1) (1998) 1-17
• [Durrleman et al. MEDIA 2008] Durrleman, S., Pennec, X., Trouvé, A., Thompson, P., Ayache, N.: Inferring brain variability from diffeomorphic deformations of currents: an integrative approach. Medical Image Analysis 12(5) (October 2008) 626-637
• [Durrleman et al. MICCAI 2008] Durrleman, S., Pennec, X., Trouvé, A., Ayache, N.: Sparse approximation of currents for statistics on curves and surfaces. In: Proc. of MICCAI'08. 390-398.
• [Durrleman et al. MEDIA 2009] Durrleman, S., Pennec, X., Trouvé, A., Ayache, N.: Statistical models of sets of curves and surfaces based on currents. Medical Image Analysis (2009).
• [Durrleman et al, MICCAI 2009] Durrleman, S., Pennec, X., Trouvé, A., Gerig, G., Ayache, N.: Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets. In Guang-Zhong Yang, David Hawkes, Daniel Rueckert, Alison Noble, and Chris Taylor, editors, Medical Image Computing and Computer-Assisted Intervention (MICCAI'09), Part I, volume 5761 of Lecture Notes in Computer Science, London, UK, pages 297-304, September 2009. Springer. Also as Research Report 6952, INRIA (06 2009)
• [Gorczowski et al. 2009] Gorczowski, K., Styner, M., Jeong, J.Y., Marron, J.S., Piven, J., Hazlett, H.C., Pizer, S.M., Gerig, G.: Statistical shape analysis of multi-object complexes. Transactions on Pattern Analysis and Machine Intelligence (2009).
• [Hazlett et al. 2005] Hazlett, H., Poe, M., Gerig, G., Smith, R., Provenzale, J., Ross, A., Gilmore, J., Piven, J.: Magnetic resonance imaging and head circumference study of brain size in autism. The Archives of General Psychiatry 62 (2005) 1366-1376
• [Khan & Beg 2008] Khan, A., Beg, M.: Representation of time-varying shapes in the large deformation diffeomorphic framework. In: Proc. of ISBI 2008. 1521-1524
• [Miller et al. 2002] Miller, M, I., Trouvé, A., Younes, L.: On the metrics and Euler-Lagrange equations of computational anatomy. Annual Review of Biomed. Eng. 4 (2002) 375-405
• [Perperidis et al. 2005] Perperidis, D., Mohiaddin, R.H., Rueckert, D.: Spatio-temporal free-form registration of cardiac MRI sequences. Medical Image Analysis 9(5) (2005) 441 - 456
• [Peyrat et al. 2008] Peyrat, J.M., Delingette, H., Sermesant, M., Pennec, X., Xu, C., Ayache, N.: Registration of 4D Time-Series of Cardiac Images with Multichannel Diffeomorphic Demons. In: Proc. MICCAI. Volume 5242 of LNCS., Springer (2008) 972-979
• [Qiu et al. 2008] Qiu, A., Younes, L., Miller, M., Csernansky, J.: Parallel transport in diffeomorphisms distinguishes the time-dependent pattern of hippocampal surface deformation due to healthy aging and the dementia of the alzheimer's type. NeuroImage 40 (2008) 68-76
• [Toga et al., 2006] Toga A.W, Thompson P.M., Sowell E.R., Mapping brain maturation, Trends in Neurosciences, 2006
• [Vaillant & Glaunès 2005] Vaillant, M., Glaunès, J.: Surface matching via currents. In: Proc. of IPMI'05. Volume 3565 of LNCS. (2005) 381-392

Image Analysis Tools

The RegistrationTool : the fast and simple image registration module of MedINRIA.

1. Introduction
2. Different Registration Solutions
3. Result Evaluations
4. Conclusions and Perspectives

Introduction

Processing and visualisation of medical data is nowadays a common yet challenging problem. The diversification and complexification of the available sources of medical information raised the need for new methods to have a better and simpler insight of these data. While clinicians need new solutions to fully exploit their data, they often cannot access state-of-the-art methods as they are mostly available in complex software targeting researchers. The MedINRIA software consists of a collection of tools that optimally exploit various types of data (e.g., 3D images, diffusion tensor fields, neural fibres as obtained in DT-MRI), providing state-of-the-art algorithms while keeping a user-friendly graphical interface. In this report we will focus on one of the applications provided by MedINRIA : the Registration Tool. This application provides a simple interface for several kinds of registration. The goal of this toolkit is to provide the clinicians with an easy way of comparing two images, either taken at different time points or from different modalities. As images are often not taken exactly in the same geometric positions and more importantly different orientations, they need to be registered in order to allow a faithful comparison. In a clinical context, registration problems that are targeted here include rigid registration for compensating the motion between acquisitions at different time points (longitudinal studies) or between different modalities, affine and deformable spatial normalisation for atlas based segmentations, and deformable registration for motion tracking in time sequences. The tool is thus divided into : first a collection of algorithms that address different registration problems, such as manual rigid registration, automatic affine registration, and deformable registration; second a collection of visualisation tools adapted to the comparison of two images, and the visualisation of deformation fields.

Different Registration Solutions Up

RegistrationTool takes two images as input, the ``fixed image'' and the ``moving image'', and offers up to now three ways of registering them:

• Manual Rigid Registration: In a few clicks, one can perform a manual rigid registration of two images. Rigid registration consists in a translation and a rotation. This is done very simply by asking the user to place in both fixed and moving images a point plus a direction. The two points are used to compute the translation, and the directions are used to calculate the angle of rotation. See this section for more details.

• Automatic Rigid or Affine Registration: We use the affine registration framework from the Insight Tool Kit http://www.itk.org to achieve this automatic matching of images. It consists in maximising the mutual information of the two images, and is embedded in a multi-resolution framework. No user input is required. This method is adapted to longitudinal studies for instance. See this section for more details.

• Automatic Diffeomorphic Demons Registration: This method provides a very up-to-date method for non-linear registration. This procedure is based on variants of the Thirion's demons algorithm for non-linear mono-modal image registration. It is fast, and ensures the invertibility of the transform, which is a very desirable feature in registration. See this section for more details.

An important note is that registration can be done iteratively: one can start with a manual rigid matching to roughly align the images, continue with an automatic affine that will feed the non-linear algorithm with a good initialisation.

For convenience, one can reset all registration steps. As software outputs, not only can the registered image be saved, but an output affine matrix file - concerning the affine and/or rigid registration - can also be saved and is compatible with the other MedINRIA modules. This matrix file can be loaded as input for the registration process in future use. On the other hand as soon as the user performs non-linear registration the output becomes a deformation field that can be saved in ITK metafile format and can also be loaded as input for future use.

Rigid Registration Up

Figure 1: This screen shot shows the two views in full screen mode where the zoom interaction has been used to isolate a specific area. A manual rigid registration is being performed (the ``Control Points'' button has been pressed). Yellow crosses are the 3D spatial points that need to correspond between images. Yellow circles define the directions (following the blue lines).
Note that the distance from the circle to the cross is not used, only the direction is taken into account.

This module provides easy-to-use manual rigid registration. The method is intuitive: the user navigates through the slices until two corresponding slices in the volume are found, the user can then choose a single point and a direction in the given slices. The RegistrationTool turns this information into a rigid matrix. The calculated rigid matrix is first applied to the moving image. Then the moving image is re-sampled to meet the fixed image size and spacing information. The process can be iterated on other slice orientations to achieve a truly 3D rigid registration.

Note that the moving image is also automatically re-sampled into the fixed image size and spacing information. The re-sampling process uses a linear interpolation (default) between pixel values but B-spine and Sinus Cardinal interpolation can be chosen.

Automatic Affine Registration Up

The user can switch between the available registration modes. By default the automatic affine mode is chosen. It consists in maximising the mutual information of the two images, and is embedded in a multi-resolution framework, based on ITK (see www.itk.org for details on the method). Three levels appear to be sufficient for a common use. No user interaction is necessary to perform this process. This operation can be executed iteratively for best results. You can also initially choose a rigid transformation as a first step (see this section) to avoid failure on large transformation, and then use the affine automatic mode to get optimal results.

Note that the output of both rigid and affine registration methods is an affine transformation matrix that can be saved for later use.

Automatic Deformable Registration Up

In this registration toolkit we also propose a solution for non-linear mono-modal image registration problems. This method provides a very up-to-date method [#!Vercauteren:ISBI:2006!#] based on variants of the Thirion's demons algorithm [#!Thirion:MIA:1998!#] to retrieve deformations between images. Although this procedure takes longer than the affine automatic mode, it allows the registration of images from different subjects, or to align an atlas image with a patient specific one. The output is a deformation vector field. We want to thank Tom Vercauteren and Mauna Kea Technologies (http://www.maunakeatech.com) for their participation on this registration algorithm.

You can save the resulting deformation field in ITK metafile format (.mha). The deformation field maps the moving image onto the fixed image. A deformation field is represented as an image whose pixel type is the 2D or 3D displacement vector. This output field can also be reload.

Image Comparison Up

Figure 2: Comparing images is easy in the RegistrationTool. This figure shows two images in comparison, the spatial positioning of both images are linked to synchronise the navigation. A full screen mode is also available (not shown here).

Comparing images remains an important goal in medical imaging. In the RegistrationTool, an effort has been made in providing intuitive and easy ways of doing so. Hence the user can simply load two images and compare them in the module. As shown in Fig. 2, spatial positioning of both images are linked together to synchronise the navigation. One can also choose a specific slice orientation and display in full screen this slice of both images.

Figure 3: This figure shows two ways of comparing images. On the left the fuse opacity option has been chosen : the two images are blended together with an alpha value the user can change. On the right the two images are turned into a grid, every second square contains alternatively the first and the second image. The size of the squares can also be changed.

As shown in figure 3, we propose two different modes. The first one, shown in Fig. 3 left, fuses the first image with the second one by blending them with an alpha value that the user can change.

The second mode shown in Fig. 3 right, uses a checker-board technique: the view is divided in squares (the number of squares can be changed), and every second square contains alternatively the first and the second image.

The navigation between slices and the change of parameters are real-time processes. Hence this is an easy method for the user to quickly compare two images, whether they are well registered to each other or not, or qualitatively evaluate their geometry differences.

Figure 4: This figure shows a comparison mode after non-linear registration. In the control panel you can choose to display the resulting deformation field (right) in parallel to the comparison view (left). The deformation field is in fact wrapped into a regular grid. Note that this feature is only available for 3D image registration.

As soon as one performs a deformable registration process (or loaded a deformation field as input), the RegistrationTool tries to visualise this field in the compare mode. The deformation field is used to wrap a regular grid. The output grid is then shown in parallel of the comparison view. Again the positions are linked to synchronise the navigation (see Fig. 4).

Conclusions and Perspectives Up

The RegistrationTool module of the MedINRIA software provides a user-friendly interface to perform simple to complex registrations, and an intuitive system for comparing two images. In Health-e-Child, the RegistrationTool is included as a service available on the HeC USB stick.