Multi-objective Optimization and Stochastic Analysis

in Focused Ultrasonic Therapy Simulation

T. Clees, N. Hornung, I. Nikitin, L. Nikitina and D. Steffes-lai

Fraunhofer Institute for Algorithms and Scientific Computing, Sankt Augustin, Germany

Keywords: Stochastic Modeling and Simulation, Multi-objective Optimization, Healthcare.

Abstract: We present new results in stochastic multi-objective optimization applied to focused ultrasonic therapy

planning. This type of non-invasive therapy uses focused ultrasound for the destruction of tumor cells and

magnetic resonance tomography for identification of tumor volume and healthy organs. During the therapy

planning the treatment parameters, such as frequency and intensity of ultrasound, are adjusted to achieve

maximal tumor destruction and minimal influence to the healthy organs. For this purpose multi-objective

optimization is used. RBF metamodeling is employed for continuous representation of discretely sampled

results of numerical simulation and for evaluation of inherent uncertainties. We apply two algorithms for

multi-objective optimization capable of non-convex Pareto front detection in the considered problem. Cross-

validation procedure and sensitivity analysis are used for estimation of uncertainties. A realistic application

case demonstrates the efficiency of the approach.

1 INTRODUCTION

Due to the non-invasive nature of the focused

ultrasonic therapy its control is often limited to

imaging methods, e.g. MRT. Numerical simulation

becomes an important step for the therapy planning.

Efficient methods for the focused ultrasonic

simulation have been presented in paper (Georgii et

al., 2011). It uses a combination of Rayleigh-

Sommerfeld integral for near field and angular

spectrum method for far field computations, which

allows determining the pressure field in

heterogeneous tissue. The bioheat transfer equation

is used to determine the temperature increase in

therapy region. Thermal dose is defined according to

CEM model or Arrhenius model (Nandlall et al.,

2009); (Pearce, 2009) as a functional of temperature-

time dependence in every spatial point in therapy

region. The simulation is considerably accelerated

by GPU based parallelization.

The purpose of therapy planning is a

maximization of thermal dose inside the target zone

(TDin) and minimization of thermal dose outside

(TDout). As usual in multi-objective optimization

(Ehrgott AND Gandibleux, 2002), the optimum is

not an isolated point but a hypersurface (Pareto

front) composed of points satisfying a tradeoff

property, i.e. none of the criteria can be improved

without simultaneous degradation of at least one

other criterion. Thus, for a two-objective problem,

the Pareto front is a curve on the plot (TDin, TDout)

bounding the region of possible solutions, see Fig.2.

Efficient methods have been developed for

determining the Pareto front.

The simplest way is to convert multi-objective

optimization to single objective one, by linearly

combining all objectives into a single target function

t(x)=

∑ w

(x) with user-defined constant weights

. Maximization of the target function gives one

point on Pareto front, while varying the weights

allows to cover the whole Pareto front. In this way

only convex Pareto fronts can be detected, because

non-convex Pareto fronts produce not maxima but

saddle points of the target function. There are

methods applicable also for non-convex Pareto

fronts.

Non-dominated set algorithm (NDSA) finds a

discrete analogue of Pareto front in a finite set of

points. For two points f and g in optimization criteria

space the first one is said to be dominated by the

second one if f

≤ g

holds for all i=1..Ncrit. A point f

belongs to non-dominated set if there does not exist

another point g dominating f. (Kung et al., 1975)

implements a fast recursive procedure to find all

non-dominated points in a given finite set.

Clees T., Hornung N., Nikitin I., Nikitina L. and Steffes-lai D..

Multi-objective Optimization and Stochastic Analysis in Focused Ultrasonic Therapy Simulation.

DOI: 10.5220/0004421900430048

In Proceedings of the 3rd International Conference on Simulation and Modeling Methodologies, Technologies and Applications (SIMULTECH-2013),

pages 43-48

ISBN: 978-989-8565-69-3

 2013 SCITEPRESS (Science and Technology Publications, Lda.)

Figure 1: MRT slices are used as a basis for volumetric material model.

Figure 2: Robust multi-objective optimization using

metamodeling of simulation results. Red points denote

simulation results, the continuous green region indicates a

metamodel of simulation results, the blue curve the Pareto

front, the black point with error bars an optimal

representative with confidence limits.

A continuous method is LP-based local

improvement algorithm (LIA). It produces a

trajectory towards the Pareto front starting from an

initial design. Each step solves a linear program

(Hornung et al., 2010)

maximize

;

where (f

,x) ≥  ≥ 0 and - ≤ x

≤ 

w.r.t. the threshold

 and the step x, for the given

gradients of objectives f

and the size of trust

region [-,]. Here i=1..Ncrit, j=1..Npar. The

approach ensures that all criteria have an

improvement at least

, maximally possible within

the given trust region. The algorithm terminates at

Pareto front, where no further improvements are

possible.

In the case when only a restricted number of

simulations is available metamodeling of simulation

results becomes useful. The objectives are

interpolated continuously in between simulation

results, and optimization algorithms are applied to

interpolated functions. In particular, RBF

metamodeling (

Buhmann, 2003) has the advantage

of generic non-degeneracy in arbitrary dimensions

and the availability of a tolerance predictor. The

interpolated function f(x) is represented as a linear

combination of special functions () depending only

on the distance to the sample points x

f(x) = ∑

k=1..Nex

(|x-x

(1)

The coefficients c

can be found from known

function values in sample points f(x

) by solving a

moderately sized linear system with a matrix 

=(|x

-x

|).

In practical problems, different sources of

uncertainty should be considered. Uncertainties of

metamodeling are related with the quality of

sampling in parameter space. These uncertainties are

reduced with the increasing density of sampling, i.e.

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

when more simulations are included into analysis.

This type of uncertainty can be estimated by a cross-

validation procedure, measuring a sensitivity of the

function value to the removal of sample points. With

the usage of Sherman-Morrison-Woodbury formula

it can be found analytically (Nikitin et al, 2012):



= f(x

) - f(x

xk remove

= c

/(

-1

)

(2)

Another type of uncertainties comes from physical

model. In addition to optimization parameters x, the

objectives also depend on a number of model

parameters p, such as physical properties of

biomaterials. These parameters are measured

imprecisely. The corresponding variation of results

can be estimated by means of the first order

reliability method FORM (Clees et al., 2012):

(f) = (∑

=1..N

(f/



(p

))

)

1/2

(3)

where non-correlated variations of parameters p

are

assumed, and () denotes standard deviation.

Finally, these two sources of scatter can be

combined into a single measure:

err = (

+

)

1/2

(4)

In further sections we apply the above described

methodology for the multi-objective optimization

and stochastic analysis of focused ultrasonic therapy

simulation.

2 MODELING

A generic workflow for modeling of the focused

ultrasonic therapy has been described in our paper

(Borsotto et al., 2012). Numerical simulation uses

FUSimlib software (http://www.simfus.de ) on 512 x

512 x 256 voxel grid. Ultrasound has been focused

in the center of the target zone for the neutral breath

state. The computational model includes angular

spectrum method with heterogeneous tissue and

reflections. CEM model is used for determination of

thermal dose. The result after 10 seconds of

exposure time has a form of spatial distributions of

pressure amplitude, temperature and thermal dose.

Fig.3 shows a typical result for thermal dose on slice

97/256 near the focal point.

The registration of tissue deformation on the

basis of MRT images is used for proper

characterization of the breathing process. It has been

shown in (Borsotto et al., 2012) that the breathing

process is one of the major sources of uncertainties

in focused ultrasonic therapy simulation. For

characterizing the breathing process two image

sequences have been processed, each containing 104

z-slices with 320 x 320 xy-resolution. One sequence

corresponds to the breath-in state, another one to the

neutral breath state. Point-to-point correspondence

between images has been determined using a block-

matching method. Fig.1 shows the resulting field of

motion vectors for slice 51/104 of the neutral breath

sequence.

Figure 3: Typical simulation result, thermal dose (range:

0eq.min/blue to 0.6eq.min/red). The target zone in the

neutral breath state is marked by the white circle.

Figure 4: Spatial material distribution: gel (black), liver

(green), other soft tissue (red), cartilage (blue), bones

(yellow). White color marks the target zone.

Segmentation procedure defines spatial distribution

of biomechanical characteristics on the voxelized

model. 5 types of materials are introduced shown

with different colors on Fig.4, for the slice 21/104. A

target zone has been modeled as a ball of 1cm

diameter located on the upper part of the liver. Its

motion according to breathing process is evaluated

on the basis of the motion vector field found in the

registration procedure and controlled by a single

parameter (t=0…1) for the transition from neutral to

breath-in state. A finer 512 x 512 x 256 voxel grid

Multi-objectiveOptimizationandStochasticAnalysisinFocusedUltrasonicTherapySimulation

has been used for simulation, where the distribution

of material has been subsampled into, using proper

positioning transformations.

3 OPTIMIZATION

Optimization is performed with respect to the

following parameters.

Table 1: Optimization parameters and their variations.

frequency 0.25…0.75 MHz

initial particle speed 0.23...0.282 m/s

The frequency of transducer is an important

parameter controlling focused ultrasonic therapy

simulation. The other one, initial particle speed, is

proportional to an acoustic intensity emitted by the

transducer (Georgii et al., 2011). As optimization

objectives the thermal dose inside and outside the

target zone have been defined as sums of the thermal

dose over corresponding voxels, ∑TDin / ∑TDout.

The variation range of optimization parameters was

regularly sampled with 25 simulations, from which

19 fall in the region of interest, shown on Fig.5. RBF

metamodel constructed on simulation results shows

Pareto front of non-convex type, for which NDSA

and LIA Pareto front detectors have been applied.

Figure 5: Results of multi-objective optimization. Red

points denote simulation results, the green points indicate

a metamodel of simulation results, the blue points indicate

the Pareto front detected by NDSA, the magenta points

show trajectories of interior points towards Pareto front

found with LIA. Black crosses show the tolerances of the

metamodel found with cross-validation procedure.

For convenience of the therapy planning it is

advisable to use a special tool for design parameter

optimization, DesParO (Clees et al., 2012), see

Fig.6. It allows to change interactively optimization

parameters and to see immediately the variation of

optimization criteria. Constraints can be set e.g.

maximizing one objective and minimizing the other,

in this way the Pareto front can be explored.

Graphical representation of interdependencies

between parameters and criteria allows to find most

influencing parameters and most sensitive criteria.

Also, the uncertainties of metamodeling found with

cross-validation procedure are shown (the red bars

under criteria sliders).

Further we focus on uncertainties coming from

physical model.

Figure 6: Optimization problem in DesParO Metamodel

Explorer. On the top: sliders of parameters, in the center:

sliders of criteria, on the bottom: a pattern of

interdependencies between parameters and criteria.

4 SENSITIVITY ANALYSIS

In our previous paper (Borsotto et al, 2012) the

sensitivity of the result to variation of 29 parameters

has been evaluated, including 5 biomechanical

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications

characteristics for 5 materials each, 3 blood

characteristics and 1 breathing parameter. It has

been shown that the parameters most influencing the

result are the absorption coefficients for gel, soft

tissue and liver as well as the breathing parameter.

For improving accuracy of the physical model a

precise measurement of absorption coefficients as

well as synchronization of transducer focal point

with breathing process has been proposed.

Here we present a new estimation of sensitivity,

based on improved measurements of critical

parameters given in (Peters, 2007).

For estimation of sensitivity, the parameters have

been varied in ± interval according to Table 2

below, then a central difference scheme has been

used to estimate entries of the Jacobian matrix:

=  crit



par

(5)

The sensitivity matrix is defined as

= J



(6)

and presented in Table 3. Then total r.m.s. of two

criteria is evaluated:

(crit

) = (∑

=1..N

)

1/2

(7)

and presented in Table 4 together with the mean

values of the criteria.

Comparing with the previous result, we see

significant improvement of precision, especially if

breathing influence is compensated. These

uncertainties should be taken into account during the

therapy planning.

Table 2: Standard deviations of absorption coefficient (in

1/m), from (Peters, 2007).

gel

0.001

liver

1.2

soft tissue

0.7

Table 3: Sensitivity matrix for objectives ∑TDin / ∑TDout

(in eq.min).

gel soft

tissue

liver

absorption -0.2

/ -0.3

-183.6 /

277.3

709 /

378.5

breath -505.2 / 546.6

Table 4: Standard deviations and mean values for ∑TDin /

∑TDout (in eq.min).

∑TDin ∑TDout

r.m.s.

total

889.5 720.4

r.m.s.

breath

compens.

732.1 469.2

mean 2273.7 3743.8

5 CONCLUSIONS

A generic approach for focused ultrasonic therapy

planning on the basis of numerical simulation has

been presented including multi-objective

optimization and stochastic analysis. Its application

to a realistic test case has been demonstrated. RBF

metamodeling of simulation results has been

performed for continuous representation of two

optimization objectives. Non-convex Pareto front of

the objectives has been determined by means of non-

dominated set and local improvement algorithms.

Uncertainties of metamodeling have been estimated

by means of cross-validation procedure. These

uncertainties can be reduced with the increasing

density of sampling, i.e. including more simulations

into analysis. Uncertainties of physical model have

been estimated by means of sensitivity analysis.

Improved accuracy of the physical model and

compensation of the influence of the breathing

process provide better precision of the result.

ACKNOWLEDGEMENTS

We are grateful to our colleagues from Fraunhofer

MEVIS: Dr. Caroline v. Dresky, Dr. Sebastian

Meier, Dr. Joachim Georgii for providing us with

MRT images, a segmented model, simulation

software as well as for fruitful discussions. This

work was supported by the Fraunhofer Internal

Program under Grant No. MAVO 821012 (FUS).

REFERENCES

Georgii, J., Dresky, C. V., Meier, S., Demedts, D.,

Schumann, C., Preusser, T., Focused Ultrasound -

Efficient GPU Simulation Methods for Therapy

Planning, in Proc. Workshop on Virtual Reality

Interaction and Physical Simulation VRIPHYS, Lyon,

France, 2011, J. Bender, K. Erleben, and E. Galin

Multi-objectiveOptimizationandStochasticAnalysisinFocusedUltrasonicTherapySimulation

(Editors), Eurographics Association 2011, pp. 119-

128.

Nandlall, S., Arora, M., Schiffter, H. A., Coussios, C.-C.:

On the Applicability of the Thermal Dose Cumulative

Equivalent Minutes Metric to the Denaturation of

Bovine Serum Albumin in a Polyacrylamide Tissue

Phantom. Proc. 8th Int. Symp. Therapeutic Ultrasound

(AIP), 1113:205-209, 2009.

Pearce, J. A.: Relationships between Arrhenius models of

thermal dose damage and the CEM 43 thermal dose.

Energy-based Treatment of Tissue and Assessment V

(Proceedings of SPIE), editor: Ryan, T.P., 7181, 2009.

Ehrgott, M., Gandibleux, X., (Eds.), Multiple criteria

optimization: state of the art annotated bibliographic

surveys, Kluwer 2002.

H. T. Kung, F. Luccio, and F. P. Preparata. On finding the

maxima of a set of vectors. Journal of the ACM,

22(4):469–476, 1975.

Hornung, N., Nikitina, L., Clees, T., Multi-objective

optimization using surrogate functions, in Proc. 2nd

International Conference on Engineering

Optimization, September 6-9, 2010, Lisbon, Portugal.

Buhmann, M. D., Radial Basis Functions: theory and

implementations, Cambridge University Press, 2003.

Borsotto D., Clees T., Nikitin, I., Nikitina, L., Steffes-lai

D., Thole C.-A., Sensitivity and robustness aspects in

focused ultrasonic therapy simulation, in J.Herskovits,

Ed., CDROM Proc. EngOpt 2012, third Int. Conf. on

Engineering Optimization, Rio de Janeiro, Brazil, July

1-5, 2012, #268.

Nikitin, I., Nikitina, L., Clees, T., Nonlinear metamodeling

of bulky data and applications in automotive design, in

M. Günther et al (eds), Progress in industrial

mathematics at ECMI 2010, Mathematics in Industry

(17), Springer, 2012, pp. 295-301.

Clees, T., Nikitin, I., Nikitina, L., Thole, C.-A., Analysis

of bulky crash simulation results: deterministic and

stochastic aspects, in N.Pina et al (Eds.): Simulation

and Modeling Methodologies, Technologies and

Applications, AISC 197, Springer 2012, pp.225-237.

Peters, K., Experimentelle Untersuchungen zur

nichtinvasiven Gewebeablation durch

hochenergetischen fokussierten Ultraschall (HIFU),

Dissertation 2007, Ludwig-Maximilians-University,

Munich (in German).

SIMULTECH2013-3rdInternationalConferenceonSimulationandModelingMethodologies,Technologiesand

Applications