Send correspondence to:
Tanja Alderliesten. Email<EMAIL_ADDRESS>
Cedric Rodriguez. E-mail<EMAIL_ADDRESS>
# Bi-objective optimization of organ properties for the simulation of
intracavitary brachytherapy applicator placement in cervical cancer

Cedric J. Rodriguez ID Dept. of Radiation Oncology, Leiden University Medical
Center (LUMC), P.O. Box 9600, 2300 RC Leiden, The Netherlands  Stephanie M. de
Boer ID Dept. of Radiation Oncology, Leiden University Medical Center (LUMC),
P.O. Box 9600, 2300 RC Leiden, The Netherlands  Peter A. N. Bosman ID
Evolutionary Intelligence Group, Centrum Wiskunde & Informatica (CWI), P.O.
Box 94079, 1090 GB Amsterdam, The Netherlands  Tanja Alderliesten ID Dept. of
Radiation Oncology, Leiden University Medical Center (LUMC), P.O. Box 9600,
2300 RC Leiden, The Netherlands

###### Abstract

Validation of deformable image registration techniques is extremely important,
but hard, especially when complex deformations or content mismatch are
involved. These complex deformations and content mismatch, for example, occur
after the placement of an applicator for brachytherapy for cervical cancer.
Virtual phantoms could enable the creation of validation data sets with ground
truth deformations that simulate the large deformations that occur between
image acquisitions. However, the quality of the multi-organ Finite Element
Method (FEM)-based simulations is dependent on the patient-specific external
forces and mechanical properties assigned to the organs. A common approach to
calibrate these simulation parameters is through optimization, finding the
parameter settings that optimize the match between the outcome of the
simulation and reality. When considering inherently simplified organ models,
we hypothesize that the optimal deformations of one organ cannot be achieved
with a single parameter setting without compromising the optimality of the
deformation of the surrounding organs. This means that there will be a trade-
off between the optimal deformations of adjacent organs, such as the vagina-
uterus and bladder. This work therefore proposes and evaluates a multi-
objective optimization approach where the trade-off between organ deformations
can be assessed after optimization. We showcase what the extent of the trade-
off looks like when bi-objectively optimizing the patient-specific mechanical
properties and external forces of the vagina-uterus and bladder for FEM-based
simulations.

###### keywords:

Bi-objective optimization, evolutionary algorithms, finite element method,
mechanical simulations, image registration, brachytherapy, cervical cancer

Figure 1: (A) Anatomical location of the vagina, cervix, uterus, and bladder.
(B) The placement and main components of an applicator, such as the tube-
shaped tandem and the lunar or capsule-shaped ovoids, used for the localized
delivery of radiation dose in BT.

##  1 INTRODUCTION

Standard treatment for women with locally advanced cervical cancer consists of
external beam radiation therapy (RT) with concurrent platinum-based
chemotherapy followed by Magnetic Resonance Imaging (MRI)-guided brachytherapy
(BT).[1] Deformable image registration could play a key role here by finding
accurate spatial mappings (tissue correspondence) between the external beam RT
treatment planning images and the BT treatment planning images of the
abdominal region. The anatomy of the relevant organs during treatment are
depicted in Figure 1. Accurate spatial mappings could potentially be used to
transfer the previously delivered dose distributions induced by the external
beam RT treatment to the BT treatment planning images, such that previous
delivered tumour and organ dosages (such as the bladder and rectum) could be
considered during the planning of the BT treatment. Finding an accurate
spatial mapping is particularly challenging in the case of cervical cancer due
to the large natural variability of the internal organs of the abdominal
region over time, such as uterine displacements and the filling of the bladder
and rectum.[2] Additionally, the intracavitary placement of an applicator that
is used for delivery of the radiation dose, as depicted in Figure 1, causes
very large deformations as well as a large content mismatch when comparing the
images acquired for BT planning to the images acquired for diagnostic purposes
and external beam RT planning.
Wide use of deformable image registration techniques in clinical practice is
still hindered due to a lack of relevant ground truth data sets. Potentially,
this could be overcome with physics-based virtual phantoms that ensure
physically plausible deformations. Many publications acknowledge the
importance of the organ mechanical properties, fixations, and external forces,
for the final physics-based deformation due to the patient-to-patient
variability. Yet, most of these works use generic mechanical properties based
on literature.[3, 4, 5, 6, 7, 8, 9, 10, 11] A potential solution is to
explicitly optimize these properties on a patient-to-patient basis. Some work
has attempted to single-objectively optimize the elastic properties of the
vagina or muscles of the pelvic floor using multi-organ Finite Element Method
(FEM) simulations of the pelvic region, showing promising results. [12, 13]
In this work, we consider finding the simulation parameter settings such that
there is an optimal deformation of both the vagina-uterus and the bladder
between diagnostic MRI scans and BT planning MRI scans. This is challenging
because these organs are directly adjacent in the abdominal region. Unless we
have the perfect simulation setup in which everything that is relevant is
modelled correctly, we hypothesize that the optimal deformation of both organs
cannot be acquired simultaneously. Hence, there might not be a single Utopian
simulation parameter setting that maximizes the deformations of both organs
that are optimal within the limits of the employed model. We therefore propose
to take a bi-objective approach where one objective is to acquire the optimal
deformation of the vagina-uterus and the other objective is to acquire the
optimal deformation of the bladder. In contrast to a single-objective
approach, a bi-objective optimization approach results in an approximation of,
not one optimal setting of simulation parameters, but a set of simulation
parameter settings (set of solutions). This optimal set is also known as the
Pareto set. This solution set, by definition, only contains non-dominated
solutions. A solution is considered non-dominated when no other solution
exists with better objective values in all the objectives simultaneously. The
objective values of each solution in the Pareto set form the Pareto front,
representing the different optimal trade-off solutions between the best
fitting deformation of the vagina-uterus and the deformation of the best
fitting bladder. This set provides insight in the range of most plausible
organ deformations, given the fidelity of the simulation. A clinical expert
can then decide which trade-off is considered most plausible or relevant.
In this work, we are the first to approach the problem of patient-specific
tuning of organ properties and external forces for simulation purposes using a
multi-objective optimization approach. Specifically, we apply an Evolutionary
Algorithm (EA) to bi-objectively optimize ten pre-selected mechanical
properties of the vagina-uterus and bladder to achieve physically plausible
organ deformations using FEM-based simulations of the placement of an
applicator that is used for delivery of the radiation dose in BT for cervical
cancer. This allows us to quantitatively and qualitatively investigate what
the extent of the trade-off between the vagina-uterus and the bladder
deformations looks like.

Figure 2: Pipeline for the optimization of the simulation parameters (i.e.,
mechanical properties and external forces of the vagina-uterus and bladder) to
improve the match between simulated and real deformations resulting from
applicator placement for brachytherapy (BT) for cervical cancer.

##  2 OPTIMIZATION PIPELINE

###  2.1 Data

We retrospectively used T2-weighted MRI scans of two patients treated with BT
for cervical cancer. All scans are acquired using the 1.5T Ingenia MRI scanner
(Philips Healthcare, Best, the Netherlands). The data of both patients
consists of an MRI scan acquired for diagnostic purposes. The data of patient
one additionally consists of two BT planning MRI scans acquired after two
separate applicator placements. During the first placement, the Venezia
gynaecological applicator (Elekta, Veenendaal, the Netherlands) was used with
an ovoid size of 26$mm$ without a cap. During the second placement, the
Utrecht gynaecological applicator (Elekta, Veenendaal, the Netherlands) was
used with an ovoid size of 20$mm$. During both placements, a tandem length of
60$mm$ at 30∘ was used. The data of patient two consists of one BT planning
MRI scan acquired after one applicator placement. Patient two was treated
using the Venezia gynaecological applicator with an ovoid size of 22$mm$ with
a cap. A tandem length of 60$mm$ at 30∘ was used. In all images, the uterus,
cervix, vagina, bladder, and abdomino-pelvic cavity were contoured for this
study by a radiation therapy technologist and checked by a radiation
oncologist. In this work, the vagina, cervix, and uterus were delineated as
one anatomical structure. We therefore denote this structure, that starts from
the vagina and ends at the uterus, as the vagina-uterus.

###  2.2 Data preprocessing

Figure 2 shows the pipeline proposed in this work. The first step in the
pipeline consists of a data preprocessing step which converts the organ
contour data into mesh representations. A model generation process converts
the contours of the organs from a 3D binary label map based on a diagnostic
MRI scan into a tetrahedron mesh. The resulting meshes are then transformed to
the desired coordinate system of the target MRI based on translations and
rotations of a rigid registration and rigid transformation process. This
process aligns the bony anatomy of the patient and is used as the initial
alignment of the diagnostic MRI scan and the BT planning MRI scan for the bi-
objective optimization. The applicator geometry is extracted from the
treatment plan and is then converted into a surface mesh. In the applicator
insertion process, we simulate the initial deformation of the vagina-uterus by
inserting the applicator mesh inside the uterine cavity while assuming generic
simulation parameters as described in Section 2.3.1. The lower bound of the
literature reported parameters ranges are used, which corresponds to low
forces and easily deformable tissue. To ensure stable bases of the organs, the
20$mm$ and 10$mm$ most inferior nodes of the vagina-uterus and bladder meshes,
respectively, are constrained (i.e., fixed) in the superior-inferior
direction. The abdomino-pelvic cavity and the applicator are considered to be
rigid and fixed objects. All the details of the data preprocessing step can be
found in Appendix A. Further, example mesh representations at different stages
are illustrated in Figure 3.

Figure 3: (A) Mesh of the vagina-uterus, bladder, and abdomino-pelvic cavity
after the model generation process but before the applicator insertion process
of placement 1 of patient 1. (B) Applicator mesh after the applicator
extraction process. (C) Mesh shape of the vagina-uterus, bladder, abdomino-
pelvic cavity, and applicator after the applicator insertion process but
before the bi-objective optimization step. The applicator mesh is contained
inside the vagina-uterus mesh. (D) Example organ mesh shape of the vagina-
uterus, bladder, abdomino-pelvic cavity, and applicator at the end of the bi-
objective optimization step.

###  2.3 Bi-objective optimization

We have chosen to use an EA as the optimization algorithm because they are
naturally suited for multi-objective optimization.[14] We specifically
consider a state-of-the-art multi-objective estimation-of-distribution
algorithm called MAMaLGaM. We provide a brief outline of its key components,
but refer the interested reader for all details to the literature.[15] Similar
to many EAs, MAMaLGaM starts with a randomly generated set (population) of
simulation parameter settings (solutions). For each solution in the
population, a complete simulation is executed to determine a solution’s
quality (called an evaluation). The resulting quality score for the bladder
and the vagina-uterus is then assigned to each solution (simulation parameter
setting). The algorithm iteratively (through generations) selects the best
solutions and generates new solutions for the population of the next
generation. Typically, the population size, maximum number of generations, and
maximum number of evaluations of the solutions are user-defined and problem-
specific. Within each generation, MAMaLGaM uses two methods in sequence to
divide the solutions into a predefined $K$ number of clusters. For the first
$D$ clusters, where $D$ is the number of problem dimensions, the best $N_{D}$
solutions within each problem dimension are grouped together. Before applying
the second clustering method, all remaining solutions are ranked based on
domination, where the first rank contains all non-dominated solutions and the
subsequent ranks are all the non-dominated solutions after removing the
lower/previous ranks. Subsequently, a fraction $\tau$ of the population (the
best-ranked solutions) is clustered into the remaining $K-D$ clusters using
balanced k-leader-means clustering in the objective space.[16] A Gaussian
distribution is estimated for each cluster and new offspring solutions are
sampled from these distributions to form the population of the next
generation. These offspring solutions can then be evaluated by running a
simulation for each new solution (simulation parameter setting). In our
application, MAMaLGaM is terminated when the predefined maximum number of
simulations has been reached.

####  2.3.1 Simulation parameter settings

To approximate the influence of the surrounding tissue on the vagina-uterus, a
constant external force is added to all the nodes of the vagina-uterus mesh,
which introduces three to be optimized simulation parameters (one value per
orthogonal direction) initialized between -1 and 1$kN$ which corresponds to a
gravitational force of 100$kg$. The resulting internal forces inside the
tissue are governed by, in our case, linear elastic equations. This relation
is governed by the elastic modulus (where a higher modulus means stiffer
material) of the vagina-uterus and is a simulation parameter initialized
between 0.001 and 100$N/mm^{2}$.[17] To ensure the fixation of the applicator
in clinical practice, the vagina is filled with pads, which results in an
expansion of the vagina. To model this effect while limiting the number of
parameters in our simulation, four sphere shaped objects are placed inside the
cavity of the vagina. The diameters of these spheres model the amount of
padding inserted in the vagina and are initialized between 0 and 40$mm$
corresponding to the size of the applicator. This results in a total of eight
simulation parameters to be optimized for the vagina-uterus mesh. The bladder
is modelled as a hollow organ with an inner and outer surface. The internal
forces are modelled using a constant elasticity between 0.05 and
11.5$N/mm^{2}$.[18] Further, the target bladder volume, which is filled using
a Foley balloon catheter before acquiring the BT planning MRI scans, is
considered to be a simulation parameter to be optimized between 0 and
600$cm^{3}$.[19]

####  2.3.2 Simulation evaluations

Simulations are executed using the FEM open-source framework SOFA (v20.12.02).
The simulation parameters are empirically defined. The simulation is divided
in 2 $\cdot$ 103 simulation steps with an internal step size of 2.5 $\cdot$
10-3 $s$ (i.e., simulating 5$s$ in simulation time) such that the simulation
is likely to be numerically stable while limiting the number of simulation
steps. The simulation consists of four meshes: the vagina-uterus mesh, the
bladder mesh, the applicator mesh, and the abdomino-pelvic cavity mesh, see
Figure 3. The latter is used to constrain the amount of anterior and lateral
movements of the bladder. Both the abdomino-pelvic cavity and applicator
meshes are rigid triangular surface meshes, since their elasticity is much
higher than the elasticity of the vagina-uterus and bladder. After each
simulation, the deformation quality (i.e., the objective values during the bi-
objective optimization step) has to be quantified for each organ separately.
The quality of the deformed organs is assessed by comparing the associated
deformed organ contours to the BT treatment planning contours using the Dice
similarity metric.[20] Some solutions correspond to simulation failure, e.g.,
due to implausible parameter combinations or fast dynamical behaviour, which
result in unstable simulations. In such cases, the worst possible Dice
similarity score of zero is returned for both organs. This is also done when
the bladder and vagina-uterus tissue are intersecting.

##  3 EXPERIMENTS

To test that the optimal deformation of both organs cannot be acquired
simultaneously, the complete pipeline is executed on three BT applicator
placements. Literature-based settings for MAMaLGaM have been used.[21] In this
case study, the number of simulations that the EA could execute is limited
(3000 evaluations) due to the long simulation time involved with each
evaluation (20 minutes). Since MAMaLGaM is a stochastic optimization
algorithm, the optimization is repeated 20 times per applicator placement to
make sure the attained results are robust. For each applicator placement
separately, the approximated Pareto sets and fronts of all runs are
aggregated. The quality of the optimization is evaluated using an independent
quality metric (root mean square). The root mean square metric indicates the
average shortest distance between the two organ surfaces and is calculated
using the open-source 3D mesh processing software MeshLab (v2020.12). To
calculate the average distance, $2\cdot 10^{5}$ points are randomly sampled on
the optimized organ mesh. Subsequently, the shortest squared distance is
determined to the target organ surface. By averaging over all these points and
taking the square root, the final metric score is determined. Finally, the
results are also assessed qualitatively by comparing the optimized organ
contours to the target organ contours. Not all deformations can be visualized
simultaneously. We therefore decided to use a subset of 5 solutions that are
distributed along the approximated Pareto front. The subset of the solutions
per placement is determined using a greedy hypervolume subset selection
method.[22] As the name suggests, this method uses the (to be maximized)
hypervolume which is a performance indicator in the case of multi-objective
optimization algorithms. Hypervolume defines the space enclosed between the
Pareto front and a user-defined reference point. The method sequentially
selects the solution from the approximated Pareto set that has the largest
contribution to the hypervolume. This typically will lead to the selection of
the solutions that mostly characterize the approximated Pareto front.

##  4 RESULTS AND DISCUSSION

Figure 4 depicts the aggregation of all approximated Pareto fronts resulting
from the twenty individual optimization runs. Higher Dice similarity scores
are found than what is observed for the initial organ shapes after the data
preprocessing step, but before the bi-objective optimization step. This shows
that there definitely is a benefit to optimizing the simulation parameter
settings. We find a trade-off front instead of one single Utopian solution for
each applicator placement. The trade-off persists, even when evaluating the
deformation quality with another metric (i.e., root mean square of the
distances between the organ surfaces).
In contrast to multi-organ optimization, single-objective approaches typically
require the selection of a similarity metric that aggregates different
optimization objectives (e.g., the average Dice similarity score of all
organs). After optimization, this will result in one simulation parameter
setting which maximizes the average Dice similarity score. This solution is
coloured in pink in Figure 4. Since the metric is defined a priori, and since
there clearly is a trade-off between the quality of deformations, this leaves
the end user uninformed about all the (potentially equally relevant) organ
deformations, and it may lead to wrong conclusions about the capacity of the
simulation. In some clinical cases, it could for example be appropriate to
prioritize the deformation quality of the vagina-uterus. Taking a bi-objective
approach makes it possible to a posteriori choose another solution along the
approximation front, i.e., one that has a higher similarity score for the
vagina-uterus at the cost of a lower similarity score for the bladder.

Figure 4: Aggregation of the approximation fronts of twenty independent runs
using the Dice similarity score of the vagina-uterus and bladder as the two
objectives during optimization of the simulation parameter settings. The
solution coloured in pink is the solution with the highest average Dice
similarity score of the vagina-uterus and bladder. Pre-optimization similarity
refers to the Dice similarity score after the data preprocessing step but
before the bi-objective optimization step. For all three placements, the pre-
optimization similarity scores of the bladder of 0.27, 0.33, and 0.22 are
outside the plotting range. The subplots depict the (independent) similarity
scores using the root mean square metric.

To visualize the deformed organs pertaining to the solutions along the
aggregated approximation front, a set of five solutions is selected, see
Figure 5, 6, and 7. In the left-right view of Figure 5 and 6, it shows that
reducing the bladder volume and simultaneously applying a force in the
inferior direction on the vagina-uterus increases similarity of the bladder
and dissimilarity of the vagina-uterus at the same time. Likewise, increasing
the bladder volume leads to improved vagina-uterus deformations but results in
an abnormal increase in the bladder in the anterior and left-right direction.
Taking a bi-objective approach makes these trade-offs explicit, enabling a
clinical expert to select which deformation is considered most plausible or
relevant. It is not surprising that the trade-off in the organ deformations
illustrated in Figure 5 and 6 are similar, since these concern optimization
results obtained for the same patient and are based on the same diagnostic
scan derived organ meshes (and their limitations). Although the trade-offs
look similar, they are not exactly the same. The difference in the highest
Dice similarity score obtained for the vagina-uterus between placement 1 and 2
is 0.035. We expect that this has to do with the unrealistic deformation of
the optimized vagina-uterus mesh near the edge of the field-of-view of the
scans. Since the vagina-uterus is undefined outside the field-of-view of the
diagnostic MRI scans, this leads to abnormal deformation around the inferior
slices of the MRI scans. These slices are clinically less relevant for the
registration in our application, since the gradients of the radiation dose
distributions are relatively small compared to those in the cervical region.

Figure 5: Orthogonal views of patient 1 (placement 1) with the target bladder
(grey), target vagina-uterus (black), and the solution contours of five
solutions selected from the approximation front. The same colour coding of
solutions is used as in Figure 4. Figure 6: Orthogonal views of patient 1
(placement 2) with the target bladder (grey), target vagina-uterus (black),
and the solution contours of five solutions selected from the approximation
front. The same colour coding of solutions is used as in Figure 4.

Figure 7 shows the trade-off between the deformation of the vagina-uterus and
the bladder of placement 1 of patient 2. In the left-right view, the
simulation parameter setting with the highest achievable bladder Dice
similarity score also results in an inferior deformation of the vagina-uterus.
Although, during this application, this inferior deformation is more extreme
compared to patient 1. When looking at the anterior-posterior view, it appears
that the vagina-uterus is used to enforce an asymmetrical deformation of the
bladder. Furthermore, the application of patient 2 has an overall lower
vagina-uterus Dice similarity score than the applications of patient 1. This
is partly due to the larger mismatch in the field-of-view between the
diagnostic MRI scan and the BT planning MRI scan. The second reason for lower
maximum achievable Dice similarity scores is due to the large mismatch in
tumour size between the MRI scans of patient 2. The organ models are based on
diagnostic MRI scans which have been acquired prior to external beam RT with
concurrent platinum-based chemotherapy. During the treatment, large tumour
regression can occur. Only applying a global force onto the vagina-uterus
cannot lead to this local shrinkage of the tumour inside the vagina-uterus.

Figure 7: Orthogonal views of patient 2 (placement 1) with the target bladder
(grey), target vagina-uterus (black), and the solution contours of five
solutions selected from the approximation front. The same colour coding of
solutions is used as in Figure 4.

Although these results are promising, several limitations of the work have to
be discussed. In this problem formulation, only two organs are considered and
no inter organ tissue is considered. Including additional abdominal organs
(i.e., bowel or rectum) or other structures (i.e., pelvis) should improve the
clinical relevance of the simulations. Moreover, future organ models should be
able to consider the tumour regression between scans or additional image
acquisition directly prior to the applicator placement could be implemented.
By subsequently increasing the number of patients, in-depth analyses would be
possible of the simulation parameter settings within (and between) different
patients. Furthermore, the field-of-view of the simulated organs should match
the field-of-view of the BT planning MRI and abnormal deformation near the
edge of the field-of-view of the simulation should be addressed. Even though
key simulation parameters settings have been selected, it is unclear how
sensitive the estimated external forces and mechanical parameter values are to
the problem formulation (i.e., organ fixations).
The bi-objective optimization step currently takes seven days due to the
computationally expensive FEM-based simulations. To evaluate one potential
setting of the simulation parameters, one full simulation is required. It
could be interesting to experiment with expensive optimization techniques
where the expensive-to-evaluate simulations are interchanged with a quick-to-
evaluate model of the simulation such that less expensive-to-evaluate
simulations are required.[23] This should result in high optimization quality
in shorter optimization time.
Finally, improved constraint handling could be included in the optimization.
The violations of some constraints, such as overlapping organs meshes, could
for instance not be penalized with the worst possible fitness values. Instead,
a degree of the violations could be considered. Furthermore, the balance
between model resolution, simulation time, and simulation numerical stability
should be considered. On the one hand, higher model resolutions lead to longer
simulation time if the same numerically stability is desired. On the other
hand, reducing the model resolution leads to shorter simulation time and
improved numerically stability, but obviously leads to reduced simulation
fidelity.

##  5 CONCLUSION

This work shows that there is a trade-off between the deformations of
inherently simplified adjacent organ models in the abdominal region in FEM-
based simulations of the placement of the applicator for BT for cervical
cancer. Understanding and acknowledging this is important, since any a priori
single-objective optimization approach may lead to wrong conclusions about the
parameters, found for the underlying biomechanical models. By taking a multi-
objective approach using an EA, we could approximate the inherent trade-off,
enabling the possibility to a posteriori evaluate which of the deformations is
considered most plausible by a clinical expert. Our results make explicit the
extent of the inherent trade-off between the similarity metrics and their
associated deformations of the vagina-uterus and bladder due to inherent
simulation model inaccuracies, while at the same time opening the possibility
to a posteriori selection of the most plausible organ deformations, thereby
increasing the probability of success in practice.

## ACKOWLEDGEMENTS

The authors thank E.C. Harderwijk (Dept. of Radiation Oncology, Leiden
University Medical Center, Leiden, The Netherlands) for his contributions to
this study regarding the data preparation. We further thank Y. Niatsetski and
R. Zinkstok (Elekta, Veenendaal, the Netherlands) for the development of a
special module to export the brachytherapy applicator geometry. The research
is part of the research programme Open Technology Programme with project
number 15586, which is financed by the Dutch Research Council (NWO), Elekta
(Elekta Solutions AB, Stockholm, Sweden), and Xomnia (Xomnia B.V., Amsterdam,
The Netherlands). Further, the work is co-funded by the public-private
partnership allowance for top consortia for knowledge and innovation (TKIs)
from the Ministry of Economic Affairs.

## References

  * [1] Cibula, D., Pötter, R., Planchamp, F., Avall-Lundqvist, E., Fischerova, D., Haie-Meder, C., Köhler, C., Landoni, F., Lax, S., and Lindegaard, J. C., “The European Society of Gynaecological Oncology/European Society for Radiotherapy and Oncology/European Society of Pathology guidelines for the management of patients with cervical cancer,” Virchows Archiv 472(6), 919–936 (2018).
  * [2] Jadon, R., Pembroke, C., Hanna, C., Palaniappan, N., Evans, M., Cleves, A., and Staffurth, J., “A systematic review of organ motion and image-guided strategies in external beam radiotherapy for cervical cancer,” Clinical Oncology 26(4), 185–196 (2014).
  * [3] Bhattarai, A. and Staat, M., “Modelling of soft connective tissues to investigate female pelvic floor dysfunctions,” Computational and Mathematical Methods in Medicine 2018, 9518076 (2018).
  * [4] Brock, K. K., Sharpe, M., Dawson, L., Kim, S., and Jaffray, D., “Accuracy of finite element model-based multi-organ deformable image registration,” Medical Physics 32(6Part1), 1647–1659 (2005).
  * [5] Chanda, A., Unnikrishnan, V., Richter, H. E., and Lockhart, M. E., “A biofidelic computational model of the female pelvic system to understand effect of bladder fill and progressive vaginal tissue stiffening due to prolapse on anterior vaginal wall,” International Journal for Numerical Methods in Biomedical Engineering 32(11), e02767 (2016).
  * [6] Courtecuisse, H., Jiang, Z., Mayeur, O., Witz, J., Lecomte-Grosbras, P., Cosson, M., Brieu, M., and Cotin, S., “Three-dimensional physics-based registration of pelvic system using 2D dynamic magnetic resonance imaging slices,” Strain 56(3), e12339 (2020).
  * [7] Diallo, M. N., Mayeur, O., Lecomte-Grosbras, P., Patrouix, L., Witz, J. F., Lesaffre, F., Rubod, C., Cosson, M., and Brieu, M., “Simulation of the mobility of the pelvic system: influence of fascia between organs,” Computer Methods in Biomechanics and Biomedical Engineering 25(10), 1073–1087 (2022).
  * [8] Dias, N., Peng, Y., Khavari, R., Nakib, N. A., Sweet, R. M., Timm, G. W., Erdman, A. G., Boone, T. B., and Zhang, Y., “Pelvic floor dynamics during high-impact athletic activities: A computational modeling study,” Clinical Biomechanics 41, 20–27 (2017).
  * [9] Peng, Y., Khavari, R., Nakib, N. A., Boone, T. B., and Zhang, Y., “Assessment of urethral support using MRI-derived computational modeling of the female pelvis,” International Urogynecology Journal 27(2), 205–212 (2016).
  * [10] Rigaud, B., Klopp, A., Vedam, S., Venkatesan, A., Taku, N., Simon, A., Haigron, P., de Crevoisier, R., Brock, K. K., and Cazoulat, G., “Deformable image registration for dose mapping between external beam radiotherapy and brachytherapy images of cervical cancer,” Physics in Medicine & Biology 64(11), 115023 (2019).
  * [11] Venugopala Rao, G., Rubod, C., Brieu, M., Bhatnagar, N., and Cosson, M., “Experiments and finite element modelling for the study of prolapse in the pelvic floor system,” Computer Methods in Biomechanics and Biomedical Engineering 13(3), 349–357 (2010).
  * [12] Silva, E., Parente, M., Brandão, S., Mascarenhas, T., and Natal Jorge, R., “Characterizing the biomechanical properties of the pubovisceralis muscle using a genetic algorithm and the finite element method,” Journal of Biomechanical Engineering 141(1), 011009 (2019).
  * [13] Mayeur, O., Witz, J.-F., Lecomte-Grosbras, P., Cosson, M., and Brieu, M., “Patient-specific simulation: non-destructive identification method for soft tissue under large strain: application to pelvic system,” in [Computational Biomechanics for Medicine ], 131–144, Springer (2019).
  * [14] Deb, K., [Multi-Objective Optimization using Evolutionary Algorithms ], Wiley (2001).
  * [15] Bosman, P. A. N., “The anticipated mean shift and cluster registration in mixture-based EDAs for multi-objective optimization,” in [Proceedings of the 12th Annual Conference on Genetic and Evolutionary Computation ], 351–358 (2010).
  * [16] Rodrigues, S., Bauer, P., and Bosman, P. A. N., “A novel population-based multi-objective CMA-ES and the impact of different constraint handling techniques,” in [Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation ], 991–998 (2014).
  * [17] Baah-Dwomoh, A., McGuire, J., Tan, T., and De Vita, R., “Mechanical properties of female reproductive organs and supporting connective tissues: a review of the current state of knowledge,” Applied Mechanics Reviews 68(6), 060801–1–060801–11 (2016).
  * [18] Martins, P. A. L., Fonseca, A. M. R. M., Santos, A., Santos, L., Mascarenhas, T., Jorge, R. M., Ferreira, A. J., et al., “Uniaxial mechanical behavior of the human female bladder,” International Urogynecology Journal 22(8), 991–995 (2011).
  * [19] Pos, F. J., Koedooder, K., Hulshof, M. C., van Tienhoven, G., and González, D. G., “Influence of bladder and rectal volume on spatial variability of a bladder tumor during radical radiotherapy,” International Journal of Radiation Oncology* Biology* Physics 55(3), 835–841 (2003).
  * [20] Dice, L. R., “Measures of the amount of ecologic association between species,” Ecology 26(3), 297–302 (1945).
  * [21] Bosman, P. A. N., “On empirical memory design, faster selection of Bayesian factorizations and parameter-free Gaussian EDAs,” in [Proceedings of the 11th Annual Conference on Genetic and Evolutionary Computation ], 389–396 (2009).
  * [22] Guerreiro, A. P., Fonseca, C. M., and Paquete, L., “Greedy hypervolume subset selection in low dimensions,” Evolutionary Computation 24(3), 521–544 (2016).
  * [23] Wang, G. G. and Shan, S., “Review of metamodeling techniques in support of engineering design optimization,” in [International Design Engineering Technical Conferences and Computers and Information in Engineering Conference ], 4255, 415–426 (2006).
  * [24] Oelke, M., Höfner, K., Jonas, U., Ubbink, D., de la Rosette, J., and Wijkstra, H., “Ultrasound measurement of detrusor wall thickness in healthy adults,” Neurourology and Urodynamics: Official Journal of the International Continence Society 25(4), 308–317 (2006).
  * [25] Panayi, D. C., Digesu, G. A., Tekkis, P., Fernando, R., and Khullar, V., “Ultrasound measurement of vaginal wall thickness: a novel and reliable technique,” International Urogynecology Journal 21(10), 1265–1270 (2010).
  * [26] Si, H. and TetGen, A., “A quality tetrahedral mesh generator and three-dimensional Delaunay triangulator,” Weierstrass Institute for Applied Analysis and Stochastic, Berlin, Germany 81 (2006).
  * [27] Powell, M. J. D., “An efficient method for finding the minimum of a function of several variables without calculating derivatives,” The Computer Journal 7(2), 155–162 (1964).
  * [28] Martin, R. S., Peters, G., and Wilkinson, J. H., “Symmetric decomposition of a positive definite matrix,” Numerische Mathematik 7(5), 362–383 (1965).

##  Appendix A Details of data preprocessing

In this appendix, the details of the model generation, rigid registration,
rigid transformation, applicator extraction, and applicator insertion are
presented.

### Model generation

The model generation process uses the clinical expert’s manually defined outer
contours of the vagina, uterus, bladder, and abdominal cavity to create the
tetrahedron meshes. Even though the vagina, uterus, and bladder consist of an
internal cavity, only the outer surface of the organs could be contoured by a
radiation therapy technologist in diagnostic MRI scans (Table 1) due to the
poor contrast between the organ cavity and the organ tissue.

Table 1: Details of the retrospective T2-weighted MRI scans where L-R denotes the Left-Right direction, A-P denotes the Anterior-Posterior direction, and S-I denotes the Superior-Inferior direction. Patient ID | Purpose scans | Number of voxels | Voxel size ($mm$)
---|---|---|---
|  | L-R | A-P | S-I | L-R | A-P | S-I
1 | Diagnostic | 288 | 288 | 24 | 0.6250 | 0.6250 | 4.0000
1 | Placement 1 | 432 | 432 | 40 | 0.5324 | 0.5324 | 4.0000
1 | Placement 2 | 432 | 432 | 40 | 0.5324 | 0.5324 | 4.0000
2 | Diagnostic | 768 | 768 | 24 | 0.3125 | 0.3125 | 4.0000
2 | Placement 1 | 432 | 432 | 35 | 0.4167 | 0.4167 | 4.4000

To be able to insert the applicator inside the vagina and uterus, requires us
to make basic assumptions about the internal cavity of these organs. The
cavity of the vagina and bladder are determined by specifying a constant wall
thickness of 1$mm$ and 2$mm$, respectively.[24, 25] The vagina wall thickness
is chosen smaller than what is reported in literature (2.7$mm$) such that the
opposing inner wall surfaces of the vagina do not overlap. The uterus cavity
is a canal which follows the centreline of the uterus starting from the
external os towards the fundus, which is created using the plug-in Vascular
Modelling Toolkit (v1.4) of the open-source software platform 3D Slicer
(v4.11). The length of the canal is chosen to be equal to the applicator
tandem length and radius, which is 60$mm$ and 1.5$mm$ for the considered
patients. The vagina and uterus are then merged into one vagina-uterus
contour. The contours are converted into triangular meshes using 3D Slicer
(v4.11).
To reduce the simulation time, the surface mesh resolutions of the vagina-
uterus, bladder, and abdominal cavity mesh are reduced (see Table 2) using the
well-established quadratic edge collapse decimation method in the open-source
3D mesh processing software MeshLab (v2020.12). The lowest resolution is
empirically chosen that still describes the internal cavity of the organs. To
ensure proper tetrahedralisation of the surface mesh, all manifold or self-
intersecting triangular surfaces are removed. Thereafter, a watertight mesh is
created by closing the obtained mesh holes. Finally, the open-source mesh
generator GMSH (v4.8.1) is used to generate a 3D tetrahedron mesh of the
hollow vagina-uterus and bladder surface meshes using the widely used Delaunay
method.[26]
The model generation process thus outputs two hollow tetrahedron meshes of the
vagina-uterus and bladder, and one surface mesh of the abdominal cavity.

Table 2: Summary of the organ and applicator mesh resolutions (i.e., number of
vertices, triangles, and tetrahedra) resulting from the diagnostic MRI scans
and BT treatment plan, respectively, of each patient generated in the model
generation process.

| Vagina | Bladder | Abdominal | Applicator | Applicator
---|---|---|---|---|---
| uterus |  | cavity | Placement 1 | Placement 2
Patient | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | 1
Vertices | 1345 | 1385 | 502 | 502 | 490 | 499 | 1083 | 923 | 1238
Triangles | 2488 | 2495 | 1000 | 1000 | 2338 | 976 | 2500 | 2500 | 2499
Tetrah. | 4726 | 4683 | 1564 | 1597 | - | - | - | - | -

### Rigid registration & model transformation

The intensity-based rigid registration process is based on the open-source
Simple ITK library (v2.1.1) which finds the initial 3D translations and
rotations of the diagnostic MRI scans with respect to the BT planning MRI
scans. Mutual information with 50 histogram bins combined with a linear
interpolator is used as evaluation metric to quantify the similarity between
the voxels intensities. The translations and rotations can then be optimized
using the Powell optimizer using Brent line search [27] with the default
settings of 100 iterations, 100 maximum line iterations, a step length of 1,
step tolerance of 1 $\cdot$ 10-6, and value tolerance of 1 $\cdot$ 10-6. This
gradient-free optimizer acquired a visually satisfactory alignment of the
pelvic bones on this set of patients, but it is possible that other gradient-
based Simple ITK methods would work just as satisfactory.
Using the acquired 3D rotations and translations, the mesh nodes of the
vagina-uterus, bladder, and abdominal cavity can now be transformed from the
diagnostic MRI coordinate system to the target BT planning MRI coordinate
system.

### Applicator extraction

A large factor influencing the deformation of the vagina-uterus is the
applicator position, orientation, and geometry. During daily clinical
practice, the applicator is reconstructed in the MRI acquired for BT treatment
planning purposes. Using an Oncentra Brachy (v4.6.0, Elekta, Veenendaal, the
Netherlands) custom plugin developed by Elekta, the applicator geometry can be
extracted in the desired position and orientation as a voxel-based contour.
Using 3D Slicer (v4.11), the applicator contour is then converted into a
surface mesh. Finally, the resolution of the surface mesh is reduced using the
quadratic edge collapse decimation method in MeshLab (v2020.12) to increase
the simulation speed.

### Applicator insertion

To combine the vagina-uterus and applicator mesh, the vagina-uterus is pre-
deformed such that the applicator fits inside the organ. To accomplish this,
an FEM-based simulation is performed with generic mechanical parameters namely
an elastic modulus of the vagina-uterus of $0.001N/mm^{2}$, elastic modulus of
the bladder of $0.05N/mm^{3}$, external forces of $0kN$ in all directions,
sphere diameters of $0mm$, and bladder volume of $0cm^{3}$. The FEM open-
source framework SOFA (v20.12.02) is utilized, where the default simulation
pipeline is used for 100 simulation seconds with a step size of 2.5 $\cdot$
10-3 $s$. These simulation resolutions ensured stable simulations capturing
the (occasionally) fast dynamic oscillations of the organ tissue while
reducing the simulation time. For the collision detection between organ
surfaces, a brute force detection approach with default contact settings is
used. The minimal proximity intersection is used where the collision between
elements is detected when they are at an alarm distance of 8$mm$ and make
contact at 0.5$mm$ distance.
The organ deformations are solved using an Euler explicit solver and a sparse
linear solver based on the LDL decomposition.[28] For the collision models, a
triangle, line, and point collision model is used with a contact elasticity
equal to the organ elasticity.
The applicator cannot directly be placed inside the vagina-uterus mesh due to
the initial rest-state of the organ. Placement of the applicator through the
vagina entry requires detailed a priori knowledge of all actions (i.e.,
translations and rotations) applied to the applicator and is not the focus of
this work. Rather, we consider finding the end state deformation of the vagina
and uterus after the placement of the applicator. Therefore, we first pre-
deform the vagina-uterus cavity in the shape of the applicator. To accomplish
this, sphere-shaped spacers are generated along the internal cavity path.
These spacers are initialized with a diameter of zero and are enlarged
linearly in simulation time until the spheres cover most of the applicator
target geometry. Simultaneously, the spacers are displaced from the initial
vagina-uterus cavity centreline to the target position of the applicator. The
spacers form a line from the tip to the base of the tandem and from the base
of the tandem to the start of the stem of the applicator, and are evenly
distributed along these two lines. After the spacers reach their target
size/position, and as a result the vagina-uterus cavity is enlarged, the
spacers are removed and the actual applicator mesh is placed inside the
vagina-uterus mesh.