# Coarse-grained simulations of DNA and RNA systems with oxDNA and oxRNA
models: Introductory tutorial

Michael Matthies  School of Molecular Sciences, Arizona State University,
Tempe, AZ 85281, USA  Matthew L. Sample  School of Molecular Sciences, Arizona
State University, Tempe, AZ 85281, USA  Fulton School for Engineering of
Matter, Transport, and Energy, Arizona State University, Tempe, AZ 85281, USA
Petr Šulc  School of Molecular Sciences, Arizona State University, Tempe, AZ
85281, USA

## ABSTRACT

We present a tutorial on setting-up the oxDNA coarse-grained model for
simulations of DNA and RNA nanotechnology. The model is a popular tool used
both by theorists and experimentalists to simulate nucleic acid systems both
in biology and nanotechnology settings. The tutorial is aimed at new users
asking "Where should I start if I want to use oxDNA". We assume no prior
background in using the model. This tutorial shows basic examples that can get
a novice user started with the model, and points the prospective user towards
additional reading and online resources depending on which aspect of the model
they are interested in pursuing.

##  1 INTRODUCTION

The fields of DNA and RNA nanotechnology use designed nucleic acid strands to
self-assemble nanoscale structures and devices. In the past forty years since
the ground-breaking work by Ned Seeman that has established the field [1],
structures of increasing complexity and sizes have been experimentally
realized. The promising applications of the field include diagnostics,
therapeutics, molecular computation and templating for nanoscale assembly [2,
3]. One of the most commonly used constructs in the DNA nanotechnology field
is DNA origami [4]. DNA origami consists of a long (7000 bases) single-
stranded DNA scaffold strand, and shorter (order of tens of nucleotides)
staple strands that connect different domains together. Computer simulations
can provide important information into the function of these devices, but also
presents multiple challenges. The substantial system sizes ( encompassing up
to tens of thousands of nucleotides) as well as the occurrence of rare events
that may be crucial to the function of the device (e.g. breaking and creation
of base pairs) present significant modeling challenges. Hence, coarse-grained
models, where each particle represents multiple atoms [5, 6], are often
employed as opposed to models with atomistic-resolution.

Among the leading coarse-grained models for DNA and RNA nanotechnology are
oxDNA and oxRNA [7, 8, 9, 10], which have been parametrized to reproduce basic
structural, thermodynamic and mechanical properties of DNA and RNA. Where
available, the models have been shown to be in good agreement with available
experimental data. They have been utilized for nanostructure design and
verification, as well as for studies of active DNA devices such as walkers, as
well as for biophysical properties of DNA and RNA [11, 5]. In this tutorial,
we describe the basics of the models, and offer a walk-through for setting up
a simulation. It is meant for users who are new to nucleic acid nanotechnology
modeling and are wondering what are the good resources to get started. It is
beyond the scope of the tutorial to cover all aspects of molecular simulation
and modeling for DNA nanotechnology and all the capabilities of the oxDNA
models and associated tools. For a more exhaustive understanding of the oxDNA
modeling tools ecosystem, we direct readers to additional references.

##  2 The oxDNA model and when to use it

OxDNA/oxRNA are coarse-grained models primarily used for simulations of DNA
and RNA nanotechnology. The models are highly coarse-grained, with 1 rigid
body representing the entire nucleotide (see Fig. 1). The beads have multiple
interaction sites, and the interaction between the sites are empirically
parameterized to reproduce DNA and RNA hairpin and duplex melting as given by
the SantaLucia and Turner nearest-neighbor models [12, 13]. The models also
capture the transition from single strands to duplex or stem formation, and
represent the mechanical properties of highly flexible single-strands, while
duplex regions behave as extensible worm-like chains, in good agreement with
known mechanical properties of DNA and RNA [7, 9, 10]. The schematic
illustration of the oxDNA model and its interactions is shown in Fig. 1. More
recently, we also introduced the ANM-oxDNA model, which represents proteins
alongside the oxDNA or oxRNA model. The description of the ANM-oxDNA model is
beyond the scope of this tutorial, and instructions on setting-up such
simulations can be found in [14, 15].

While detailed description of the model’s parametrization and their properties
of provided elsewhere [7, 10, 9, 8], the novice user should be aware of the
following: 1) The oxDNA model comes in two versions, oxDNA1 and oxDNA2. The
latter, oxDNA2, is the recommended choice as it includes salt effects using
the Debye-Huckel potential and accounts for major and minor grooving, which
oxDNA1 did not capture. Similarly, the oxRNA2 also includes the Debye-Huckel
potential as opposed to oxRNA1, which did not support it, and hence is also
the recommended interaction to use. 2) Both models can be used with two
options: sequence-independent and sequence-dependent. The sequence independent
version is the default version. While only A-T (A-U) and C-G base pairs can
form in the sequence-independent model, their strength is the same: The model
uses averaged value for all of them. This model is suitable for studies where
one does not want the results to be occluded due to a particular choice of
interaction strength, or where the sequence effects average out due to similar
frequency of use of AT and GC base pairs. The sequence dependent model uses
stacking and base pairing interaction strengths parameterized to the melting
temperature predictions of duplexes. Additionally, in the oxRNA model, it also
allows formation of wobble (G-U) base pairs. 3) OxDNA is the name of the
coarse-grained DNA model, but is also used interchangeably to refer to the
oxDNA simulation package [16], which is a molecular simulation software
package that currently implements over 20 different coarse-grained models,
including oxDNA and oxRNA. Originally, the simulation package only supported
the oxDNA model, hence the names were synonymous. However, since the package
was first introduced as an open-source software in 2012, new models have been
added to it. Additionally, the oxDNA and oxRNA models have also been
implemented in the general-purpose LAMMPS simulation package [17]. In this
tutorial, we only focus on the simulation of oxDNA with the oxDNA simulation
package, which also includes additional set of scripts in Python for data
evaluation.

Figure 1: A schematic overview of oxDNA and oxRNA models. a) A superposition
of coarse-grained and atomistic nucleotide representations. Even though for
visualization purposes we draw the hydrogen bonding and stacking base site as
an ellipsoid, the model represents them as spherical excluded volume sites. b)
A schematic of a DNA duplex, along with a representation of a single base with
position vectors $\mathbf{b}$ and $\mathbf{n}$, as well as a schematic
illustration of the interactions in the model (Debye-Huckel repulsion between
backbone sites is not shown). c) A schematic of a melting transition in the
oxRNA model.

Prior to starting any simulations,it is important to consider if the oxDNA
tool can provide valuable insight into the problem you are studying. Namely,
one should ask if the same research problem can be answered with a simpler and
faster method, and if the process one is interested in is accurately captured
by the model, or would answering it require a more detailed / accurate
description of the system? Often, one can use secondary structure prediction
tools such as ViennaRNA or NUPACK [18, 19] for fast enumeration of possible
interaction states. As opposed to these tools, oxDNA provides a three-
dimensional representation, can be used to study kinetics, and simulation of
complex 3D structures including pseudoknots. On the other hand, the model
lacks atomistic-level resolution and cannot represent interactions with other
biomolecules, which limits its applicability to studies of complex biological
environments.

Currently, OxDNA and oxRNA are the only models for DNA and RNA nanotechnology
that can correctly and efficiently capture single-stranded to double-stranded
transitions. This makes them especially useful to study systems where base-
pair creation or breaking can occur, as well as for structures which consist
of both single-stranded and double-stranded regions. The model has been shown
to be able to study systems consisting of up to 1 million nucleotides [20].
However, simulations of diffusion-limited processes that rely on the diffusion
of single-stranded DNA or RNA, can be quite time-consuming. Currently, the
model is not fast enough to e.g. simulate assembly of DNA origami structure
from individual strands, and advanced biasing methods and significant
computational resources need to be applied to study systems of such sizes
[21].

The simulation package implements both CPU and GPU molecular dynamics engines,
and CPU implementations of the Monte Carlo and cluster Virtual Move Monte
Carlo (VMMC) algorithms. The latest version of the code is freely available at
https://github.com/lorenzo-rovigatti/oxDNA. The online repository features
code documentation, which includes installation instructions and also provides
a more comprehensive explanation of the model and its parameters. In the
following, we walk through the steps of setting up and evaluating an example
simulation. It is beyond the scope of this tutorial to cover all aspects of
molecular simulations, and for the general techniques and explanation of
theory behind molecular simulation, we recommend popular textbooks such as
[22, 23]. A more condensed introduction to molecular simulation algorithms is
also provided in [24].

##  3 Setting-up and running simulations

###  3.1 Creating the starting configuration

The first step in setting-up a simulation is the preparation of a system in
the oxDNA model file format. The oxDNA (or oxRNA) models requires two files:
the topology file (typically ending with .top suffix) and configuration file
(typically saved as a .dat file). Below, we list an example of the topology
file for a system consisting of two DNA strands, each of length 3 bases:

    6 2
    1 T -1 1
    1 A 0 2
    1 C 1 -1
    2 G -1 4
    2 T 3 5
    2 A 4 -1

where the first line indicates there are 6 nucleotides in two different
strands. Then, the sequence is listed in the $3^{\prime}$ to $5^{\prime}$
order, with one nucleotide per line. The first column is the id of the strand
the nucleotide belongs to, then the base type, and the id of the $3^{\prime}$
and $5^{\prime}$ neighbor. Alternatively, it is possible to replace the base
type (A,T,C,G) by a number with absolute value greater than 10, in which case
it can only bind to a base with id such that their sum is 3. For example, in
the following topology the first nucleotide can be only able to bind to the
last one, and no other base pairs are possible for it, thus avoiding some
possible undesired mismatches that you want to prevent from forming:

    6 2
    1 13 -1 1
    1 A 0 2
    1 C 1 -1
    2 G -1 4
    2 T 3 5
    2 -10 4 -1

Future stable release of the oxDNA tool will also support a new topology
format where the sequence is simply listed as one strand per line, but the
above format will remain supported for backward compatibility. An example of
the first three lines of the configuration file are

    t = 0
    b = 20.0 20.0 20.0
    E = -1.27  -1.42 0.14

where the first line is the timestep of the simulation (starting at 0), the
second line is the dimensions of the simulation box edges, and the third line
is the total, potential, and kinetic energy per particle respectively. It is
advised that the box size is large enough such that the simulated structure
cannot interact with itself (e.g. has to be larger than the largest distance
between any pairs of nucleotides in the system). This can be easily verified
by visualizing the structure in oxView tool [14] and switching on the Box
option in View tab. Then, after the top three header lines, the configuration
file has one line per particle (not shown above). Each line lists, for each
nucleotide in separate columns, x, y and z position, the three components of
the particle vectors $\mathbf{b}$ (pointing from the center of mass of
nucleotide to base site) and $\mathbf{n}$ (perpendicular to $\mathbf{b}$ and
will point towards $5^{\prime}$ neighbor in a formed duplex, see Fig. 1b).
Then, the three components of the velocity and angular velocity vector are
listed as well. Once the simulation is run, the simulation code will also
produce a trajectory file, which is a series of configuration files (in the
same format as described above) concatenated one after another.

It would be too tedious to manipulate and create the topology and
configuration files directly, so a set of tools are provided to convert from
common design tools such as CADNano, Tiamat or PDB format [25, 26]. Some other
existing tools, such as MagicDNA or sCADNAno [27, 28], also allow directly
saving designed system in oxDNA format.

The standard tool for oxDNA visualization, oxView, is available as a browser-
based application at https://oxView.org. It can import from common DNA
nanotechnology formats directly into oxDNA, and then save the system as .top
and .dat files in the format above. Importing structures from other design
tools is advised for origami sized systems, while smaller systems (like tiles)
can be created directly using the oxView editing interface, which allows
creation and deletions of single-stranded and double-stranded sections and
connecting them together. The oxView tool along with its use for structure
import and generation is described in detail in [14]. The conversion tool
between different formats is also provided as a standalone website
http://tacoxdna.sissa.it/ [29].

###  3.2 Installation of the oxDNA tool

The software tool is intended to be run in a high performance computing (HPC)
environment, such as a dedicated workstation or computer cluster. Such systems
typically run using Linux systems, and as such, our tool is currently intended
to be run in a Linux environment. It can be also run in a Linux virtual system
on Windows. To compile the executable (binary) of the program, go to the oxDNA
directory after downloading the tool and type

        mkdir build
        cd build
        cmake ../ -DCUDA=1
        make

If successful, the above command will lead to a binary file called oxDNA in
the build directory. The above commands also assumes you have the CUDA toolkit
library installed on your system that supports GPU computing. If you leave out
the -DCUDA=1 parameter, the code will only be compiled with support for a
single CPU simulations. In most HPC settings, one needs to import a NVIDIA
CUDA compiler and a compatible C++ compiler (we recommend GNU C Compiler)
using the module add command. Please refer to your computer cluster
documentation for details how to import the right libraries.

Optionally, you can install Python binding libraries as well, which gives you
access to the oxpy modules for preparing and running the simulations from a
Python interface, as well as providing access to the oxDNA Analysis Tools
modules for analysis of simulations in Python. In that case, the installation
process is

        mkdir build
        cd build
        cmake ../ -DCUDA=1 -DPython=1 -DOxpySystemInstall=1
        make
        make install

The above Python binding will only work with Python version 3.9 or higher. For
future reference of troubleshooting of known installation issues, please refer
to the online documentation of the tool. For additional information and
troubleshooting, follow the instructions on the github repository
documentation page https://lorenzo-rovigatti.github.io/oxDNA/.

###  3.3 Use of online resources to run oxDNA

As an alternative to installing and compiling your own oxDNA code locally,
online tools are available to perform simulations for you. To simplify the
setup and use of the model, we maintain a free online webserver
https://oxDNA.org, where the topology and configuration files can be directly
uploaded. The server is equipped with 8 GPU cards and supports oxDNA, oxRNA
and ANM-oxDNA simulations. It is meant for users that do not have access to
their own computer clusters and are primarily interested in running unbiased
simulations of DNA or RNA systems. The details of how to use this service is
provided in [30]. For certain applications this might be a better option, as
the server interface is optimized for typical use cases of the oxDNA and oxRNA
models. It provides basic online analysis tools that can produce mean
structures, preform flexibility analysis, create videos of simulated
trajectories, show bond occupancy, and calculate distributions of angles
between different helices and distances between groups of nucleotides. In the
rest of the tutorial, we assume that you are running all simulations on a
Linux cluster, and hence need to setup all simulations by yourself.

###  3.4 Setting up simulation parameters

####  3.4.1 Molecular Dynamics on CPU

Once you have prepared the starting configuration and topology, we need to
setup a text file with instructions to run the simulation. This file is
typically named input, and contains parameters of the simulations. An example
of such a file is provided below:

    ####  PROGRAM PARAMETERS  ####
    backend = CPU
    seed = 4982

    ####    MD PARAMETERS    ####
    sim_type = MD
    newtonian_steps = 103
    diff_coeff = 2.50
    thermostat = john
    dt = 0.003
    verlet_skin = 0.05

    ##### MODEL PARAMETERS #####
    interaction_type = DNA2
    use_average_seq = no
    seq_dep_file = oxDNA2_sequence_dependent_parameters.txt
    salt_concentration = 1.0
    T = 25C

    ####    INPUT / OUTPUT    ####
    steps = 1e8
    topology = 75deg.top
    conf_file = 75deg.dat
    trajectory_file = trajectory.dat
    refresh_vel = 1
    time_scale = linear
    restart_step_counter = 1
    energy_file = energy.dat
    print_conf_interval = 1e5
    print_energy_every = 1e4

Any line that starts with # are treated as a comment and are ignored by the
program. The simulation parameters are grouped by their role. First, backend
specifies we will run the code on CPU, and seed initializes the random number
generator. Omitting the seed parameter prompts the program to default to a
randomly generated number when launched. If you run one structure in two
different simulations with the same seed, they will produce identical
trajectories, so make sure to launch simulations with unique seeds. Next, we
set the parameters of the molecular dynamics (MD) algorithm. While detailed
description of molecular dynamics is available elsewhere [22], it can be
approximately understood as numerical integration of Newton’s equations of
motion. This involves parameters like timestep dt and verlet list, which is
important for neighbor list calculation, as elaborated in [22]). Each
nucleotides movement in the simulation is dictated by forces equal to the
derivatives of the oxDNA model potential. Additionally, it is coupled with a
thermostat, which ensures that the states sampled by the system correspond to
the specified temperature. In the example above, we use an Andersen-like
thermostat (specified by parameter john in reference to the article which
introduced it [31]). This thermostat periodically resamples random velocities
and angular momenta of nucleotides so that they correspond to the Boltzmann
distribution at the specified temperature. The parameters of the thermostat
are diff_coeff which sets the diffusion constant of a single nucleotide, and
newtonian_steps which is related to the frequency with which the velocities of
nucleotides are refreshed. It is advised that these parameters are kept to
their values specified above. The steps parameters specifies how long the
simulation should run (i.e. how many steps it will perform). This number is
system-dependent, and it should be as large as necessary to sample all
relevant states of the simulated system.

The model parameters section specifies that we use the oxDNA2 model version
(interaction_type = DNA2), with its sequence-dependent variant. The file
oxDNA2_sequence_dependent_parameters.txt must either be copied from the main
oxDNA directory to the directory where the simulation is launched or the full
path to the file location can be provided. The above configuration sets the
sodium concentration to 1M using the Debye-Hückel potential, and the
temperature to 25 Celsius degrees. An identical set of parameters would also
apply to the oxRNA model, however the interaction is called RNA2.
Additionally, there is an optional (but highly recommended) parameter
mismatch_repulsion = 1, which introduces a weak repulsion between mismatched
base pairs to correct for the high stability of mismtached base pairs in the
oxRNA model [9].

Finally, the input/output section specifies the name of the topology and
starting configuration files (75deg.top and 75deg.dat in this example). The
generated trajectory will be saved to trajectory.dat (specified to be saved
every $10^{5}$ simulation steps). The energy per particle as a function of
time will be saved to a file named energy.dat every $10^{4}$ steps, using
linear time scale. Setting restart_step_counter to 1 means that every time the
simulation is stated, it deletes the energy.dat and trajectory.dat files and
then restarts the simulation from time 0. If the parameter is set to 0, the
simulation will instead continue appending to existing trajectory.dat and
energy.dat files, resuming from the time specified in the configuration file.
The simulation saves the last configuration when it ends (or the job is
cancelled or killed by the server) to a file named last_conf.dat. If you use
the option restart_step_counter = 0, it is vital to ensure you are restarting
the simulation from the last simulation state by setting conf_file =
last_conf.dat.

To run the simulation, it is advised to create a separate directory which
contains the topology, configuration and input files. The simulation can then
be launched by

    $OXDNAPATH/build/bin/oxDNA input

where $OXDNAPATH is set to the oxDNA installation directory. The program will
then run and produce the output files in the same directory. To ensure you
simulation is running as intended, it is advised to observe the values of
energy printed on the screen (ideally they should be about $-1.4$ per
nucleotide for a system consisting of mostly base-paired nucleotides), and to
use oxView to visualize the produced trajectory.

####  3.4.2 Molecular Dynamics on GPU

Single CPU simulations are usually too slow for most systems of interest. If
the size of the studied system exceeds about 300 nucleotides, the GPU
simulation can provide up to two orders of magnitude speed-up. To use a GPU-
version of oxDNA, the first two sections of the input file have to be changed
to

    ####  PROGRAM PARAMETERS  ####
    backend = CUDA
    backend_precision = mixed
    use_edge = 1
    edge_n_forces = 1

    ####    MD PARAMETERS    ####
    sim_type = MD
    CUDA_list = verlet
    verlet_skin = 0.2
    max_density_multiplier = 10
    CUDA_sort_every = 0
    newtonian_steps = 103
    diff_coeff = 2.50
    thermostat = john
    dt = 0.003

The detailed explanation and fine-tuning of the parameters above are given in
the oxDNA code documentation. For a majority of large systems (e.g. DNA
origami), these listed values should provide optimal speed on the latest GPU
cards. The remaining sections setting-up the model and input and output files
are the same as they were for the CPU implementation.

####  3.4.3 Monte Carlo Simulation on CPU

The simulation code also supports both the Metropolis Monte Carlo Algorithm
and Virtual Move Monte Carlo (VMMC) algorithm for the oxDNA and oxRNA models.
The Monte Carlo algorithm relies on proposing a trial move, either with a
single particle or a group of particles, and subsequently accepting or
rejecting it [22]. The Virtual Move Monte Carlo algorithm [32] is particularly
efficient at sampling the states of smaller systems (ideally less than 100
nucleotides), which it achieves faster than molecular dynamics. The input file
for the MC algorithm (which replaces the MD PARAMETERS section of the CPU
version of the code) is:

    ####    MC PARAMETERS    ####
    sim_type = MC
    ensemble = NVT
    delta_translation = 0.10
    delta_rotation = 0.25
    list_type = cells

The parameters delta_translation and delta_rotation set the mean size of
proposed translation or rotation moves, and should ideally be chosen so that
about 30% of suggested moves are accepted. Large scale moves will be rejected
too often, while short moves will be accepted, but overall it will take too
many proposed moves before the system significantly changes. Provided
list_type = cells is one of the possible means (along with verlet lists) that
algorithm can efficiently keep track of its neighbors. For the VMMC
simulation, the recommended input parameters are

    ####    VMMC PARAMETERS    ####
    sim_type = VMMC
    ensemble = NVT
    delta_translation = 0.10
    delta_rotation = 0.25
    verlet_skin = 1.0
    small_system = 1
    maxclust = 30

In the above, there are two optional parameters: maxclust and small_system,
which in certain scenarios (typically for system sizes below 50 nucleotides)
can increase the speed of the algorithm. When the simulation tool stops, it
prints out benchmarking information (such as CPU time per one simulation step)
which can be used to fine-tune parameters for maximum efficacy.

###  3.5 Relaxation of a structure

When using designs exported directly from design tools, it is often necessary
to perform a "relaxation" simulation prior to the production run (the
simulation intended to collect data about the system). This relaxation phase
is crucial because design tools typically place nucleotides into positions
that are not physically possible. For example, adjacent nucleotides on the
same strand might be placed too far from each other, exceeding the permissible
distance that the covalent bond between nucleotide backbones would allow. The
oxDNA model force-field would in such a case result in a nearly infinite force
trying to pull these "stretched" nucleotides together, causing numerical
instability and program error. Another common issue with imported designs is
that sometimes nucleotides are placed too close to each other, resulting in
them overlapping. Since the repulsion forces in the oxDNA model will be
enormous if two nucleotides are too close to each other, this can result in
very high forces and subsequent numerical instability and program crash. It is
hence recommended to perform relaxation simulations prior to production runs.
We will list the parameters needed to run a relaxation simulation below (they
need to be added to either CUDA or CPU simulation input file):

    ####  RELAX PARAMETERS  ####
    dt = 0.001
    max_backbone_force = 10

Where dt is set to a smaller value as compared to the production run to
increase numerical stability and max_backbone_force is the parameter
permitting overstretched backbone bonds. The relaxation steps can be run for a
certain number of steps (typically $10^{5}$ to $10^{7}$, but can be longer for
larger systems with many stretched regions). In some cases (which feature lots
of overlapping nucleotides), it might also be necessary to first run Monte
Carlo simulation before proceeding to molecular dynamics relaxation. The
automated relaxation is also implemented in the oxDNA.org webserver.

###  3.6 Running simulations using the Jupyter notebook interface

Beyond the command line interface, oxDNA also supports python bindings. This
capability facilitates both oxDNA simulation and analysis directly from a
Jupyter notebook. The setup is described in more detail in the online
documentation at https://lorenzo-rovigatti.github.io/oxDNA/oxpy/index.html#an-
example-of-a-simple-simulation. To streamline the setup of routine
simulations, we provide basic boiler plate code. The Jupyter notebook code
provided below preforms the setup of a default production run simulation and
provides a convenient way to view and analyse the trajectory.

⬇

# visualisation library

from oxDNA_analysis_tools.UTILS.oxview import from_path

# boilerplate code to setup and monitor simulations

from oxDNA_analysis_tools.UTILS.boilerplate import setup_simulation,\
get_default_input,\ Simulation

\par\par\par# get default settings for a CUDA simulation

parameters = get_default_input()

# setup sequence dependence

parameters["use_average_seq"] = "no"

parameters["seq_dep_file"] = "../../oxDNA2_sequence_dependent_parameters.txt"

\par# setup simulation

# we will simulate the 75 degree layered tile from

# https://doi.org/10.1021/jacs.8b07180

input_file = setup_simulation("./75deg.top", # topology

"./75deg.dat", # starting configuraiton

"./simulation", # output path

parameters, # simulation parameters

kill_out_dir = True)

\par# create a simulation class

s = Simulation(input_file)

\par# view initial configuration

s.view_init()

\par# run the simulation

s.run()

This code will setup and run a production run of the 75 degree layered
crossover tile. To monitor the progress of the simulation you can utilize the
following code.

⬇

# visualize simulation progress

s.plot_energy()

This will produce a plot of the total energy per simulation step. The red line
indicates the completion of the simulation. This plotting function can be
rerun as many times as needed to follow the simulation progress. To visualize
progress one can additionally utilize following code, visualizing the last
configuration:

⬇

# displays the last configuration

s.view_last()

Finally after the simulation finishes one can analyze the resulting trajectory
using the following jupyter code, calling the oxDNA analysis tools mean
script.

⬇

# after the simulation is finished,

# lets perform a mean structure analysis using oxDNA_analysis_tools

oat mean ./simulation/trajectory.dat -o ./simulation/mean.dat -d\
./simulation/devs.json

The resulting files can be visualized in the Jupyter notebook using:

⬇

from_path("./simulation/75deg.top",

"./simulation/mean.dat",

"./simulation/devs.json")

We provide this code and the configuration files as a Jupyter notebook in the
example folder "BOILERPLATE_Jupyter", shipped with the oxDNA code.

##  4 Advanced simulations

The examples above covered the basic setup for unbiased simulations on GPU and
CPU. For many systems, this might not be enough: If, even after few days of
simulations, you have still not seen the simulated system sample the relevant
states often enough, this suggests that the unbiased simulation is not
providing sufficient data to be representative of the full behavior of the
structure. For example, from experimental FRET measurements, you know that the
arms of a DNA nanostructure sample certain ranges of possible distances.
However, if the oxDNA simulations do not show the entire range of the motion,
it means that it has not been run long enough, or that more advanced methods
to enhance the sampling need to be employed. The oxDNA simulation implements
discrete umbrella sampling in combination with the VMMC algorithm, which is an
ideal tool to sample the entire free-energy landscape as a function of numbers
of base pairs in the system. Examples of use include simulations of melting of
duplexes, hairpins and pseudoknots and strand displacement system [33, 34].
This method is documented, along with provided example codes, on the github
repository of the oxDNA code. We also support umbrella sampling as a function
of continuous parameter (typically distance between groups of nucleotides)
that can be used to enhance the sampling of large scale conformational change
of different structures. The example code to setup such simulations is
available as part of [35] and is also discussed briefly below. Finally, the
model also support enhanced kinetics sampling to extract e.g. rates of a
particular reaction such as strand displacement. The forward flux sampling
method is supported both on CPU and GPU in the oxDNA code, with examples
provided in the the oxDNA documentation on github.

###  4.1 Continuous parameter umbrella sampling

An open-source Python wrapper of the oxDNA Python bindings, available at
https://https://github.com/mlsample/ipy_oxDNA, provides a semi-automated
Jupyter notebook interface to run continuous parameter umbrella sampling. This
runs on the GPU implementation of the oxDNA molecular dynamics engine, and as
such allows for enhanced sampling of origami-sized systems or larger. The
Jupyter notebook interface enables a greater ease of use and supports
reproducible and open source simulation notebooks. Furthermore, NVIDIA’s
multiprocessing service can be harnessed to increase simulation throughput
over 2.5x by allowing over 40 origami sized simulations to be run on a single
GPU.

Umbrella Sampling works by adding an external biasing potential to the system
to enforce sampling all values of the chosen order parameter. An order
parameter is any measurable value that is of interest in a simulation. An
example of a chosen order parameter can be the distance between two groups of
nucleotides on the ends of a DNA origami that we expect to undergo large scale
conformational change over a range of distances in the simulation. The biasing
potential effectively ’flattens’ the energy landscape, making it easier for
the system to overcome energy barriers and explore a wider range of states. To
ensure we span the entire range of your order parameter of interest, we run
multiple simulations (windows), where each window samples only a specific
subset of the order parameter.

The set of biased simulations obtained by running different windows are then
combined using the Weighted Histogram Analysis Method (WHAM) [36]. WHAM
effectively ’removes’ the bias imposed by the external biasing potential
applied during the simulations, and then stitches the distributions of each
window together to give a single unbiased estimate of the free energy profile
for your order parameter. Careful selection of the order parameter is crucial,
as it should be able to effectively distinguish between the different states
or conformations of the system.

Following the import of the necessary python modules:

⬇

# Import necessary modules

from umbrella_sampling import ComUmbrellaSampling

from oxdna_simulation import SimulationManager

the first step to use the python wrapper is to specify the directory that
contains the oxDNA conformation and topology files (file_dir), which is
ideally prepared in the command line prior to running the simulation. Next,
you need to specify a unique name for the directory which will contain our
umbrella sampling system (system). We can then define our order parameter,
where in this example we use the distance between the center-of-mass between
nucleotides in group one (com) and group two (ref), defined by their
nucleotide index.

⬇

# Define the directory where your files are stored

file_dir = ’abs/path/to/files’

\par# Define the name for your system

system = ’user_defined_system_name’

\par# Define the particles for center-of-mass and reference

com_list = ’1,2,3’ # This is a comma-separated list of particle indices

ref_list = ’4,5,6’ # These indices can be found using oxView

We next set the umbrella simulation parameters which are determined a priori
or refined over umbrella simulation trials. These include the order parameter
sampling range (xmin, xmax), the number of simulation windows (n_windows), and
the stiffness of the external biasing potential (stiff).

⬇

# Define the parameters for the umbrella sampling

xmin = 0 #Lowest com distance attemped sampling

xmax = 10 #Highest com distance attemped sampling

n_windows = 20 #Parallelizable simulations sampling subsets of range

stiff = 1 #Stiffness or ’k’ value of umbrella potential

For the oxDNA simulation parameters, we need to set the number of steps for
each window, the temperature, and the intervals at which to print the energy
and configuration data. We can also modify other oxDNA parameters as required,
such as salt concentration.

⬇

# Define the parameters for the equilibration phase

equilibration_parameters = {’steps’:’1e7’, ’T’:’20C’,

’print_energy_every’: ’1e6’, ’print_conf_interval’:’2e6’}

# Define the parameters for the production phase

production_parameters = {’steps’:’2e7’, ’T’:’20C’,

’print_energy_every’: ’2e6’,’print_conf_interval’:’1e7’}

\par# Define a list of shared parameters

shared_params = [simulation_manager, n_windows, com_list,

ref_list, stiff, xmin, xmax] #optional list to reduce reduencay

We create an umbrella sampling object and a simulation manager object. If a
system directory doesn’t exist, instantiating a ComUmbrellaSampling object
will create one. If a system directory already exists, it will read the
current state of the directory, and will "reattach" to the directory. This
means if we need to clear the Jupyter kernel or start a new session with a
umbrella simulation run previously, we can "reattach" to the umbrella system
directory and preform analysis or continue simulations without modifying the
contents or overwriting files. The simulation manager allocates the simulation
to the available GPU and CPU resources optimally within on a single node.

⬇

# Initialize the umbrella sampling and simulation manager objects

us = ComUmbrellaSampling(file_dir, system)

simulation_manager = SimulationManager()

Before the simulations are built, we start the NVIDIA multiprocessing service
for better performance.

⬇

# Enable the Nvidia multiprocessing service

simulation_manager.start_nvidia_cuda_mps_control()

Now, you’ll prepare and run equilibration simulations. These will generate
your simulation windows to their assigned sampling range. It is recommended to
start with a pre-equilibration check, a brief 50,000 step simulation, to
verify your system setup.

If the pre-check is successful, re-run the ’build_equilibration_runs()’
method. This will overwrite your previous equilibration setup, emphasis on
rebuilding will delete previous files, and start the true equilibration
process. The number of steps for this phase depends on your system and order
parameter, but for origami-sized systems, a minimum of 2 million steps is
suggested.

The ’simulation_manager.run()’ method executes these equilibration simulations
as a set of ’grandchild’ processes. This lets you continue using your notebook
while simulations are running. However, if you prefer to block your kernel
while running umbrella simulations (necessary for a umbrella sampling python
script with Jupyter-validated umbrella setups), use the
’simulation_manager.worker_manager()’ method instead. This spawns the
simulations as a set of ’child’ processes.

⬇

# Build and run the equilibration simulations

us.build_equlibration_runs(*shared_params, equilibration_parameters)

simulation_manger.run() # This runs the equilibration simulations

As the simulations run (or afterwards), you can analyze the time series order
parameter values using the following commands:

⬇

# Analyze the simulation results while they’re running or afterwards

us.com_distance_observable(com_list, ref_list)

plt.figure(dpi=200)

for idx in range(0,n_windows,1):

us.analysis.view_observable(’eq’, idx, sliding_window=False)

plt.legend(fontsize=4)

You can view the final conformation of any window by running the following:

⬇

# View the last conformation of any window

us.equilibration_sims[idx].analysis.view_last()

Finally, you can proceed to run the production simulation and the Weighted
Histogram Analysis Method (WHAM) analysis:

⬇

# Build and run the production simulation and WHAM analysis

us.build_production_runs(*shared_params, production_parameters)

simulation_manger.run()

\parwham_dir = os.path.abspath(’../wham/wham’)

n_bins = ’200’

tol = ’1e-5’

n_boot = ’30’

us.wham_run(wham_dir, xmin, xmax, stiff, n_bins, tol, n_boot)

Finally, after you have run the WHAM technique, you can view your free energy
profile by running:

⬇

# View your free energy profile

us.plot_free()

A plethora of Jupyter notebooks and python scripts are provided in the github
repository as well.

##  5 SIMULATION ANALYSIS

After simulations have been completed, it is often necessary to extract the
relevant information using statistical analysis. Each case is system-specific,
but we provide towards this goal oxDNA Analysis Tools (oat), a Python-based
suite designed for analyzing oxDNA simulations. These tools cover some typical
use-cases for analysis of simulations, and can also be modified by the users
to suite their current needs. They can be used as a series of standalone
command-line scripts or as importable Python modules for creating personalized
analysis tools. The command-line scripts can be accessed by the command oat
<script name> <script arguments>, with autocomplete features available for
script names. The suite includes a diverse range of scripts for various
functions like bond analysis, clustering, and principal component analysis,
among others. Additionally, the scripting API for oat can be imported into
Python scripts, providing utilities such as the optimized oxDNA trajectory
reader. Comprehensive API documentation is available for both top-level
analysis and utility APIs at https://lorenzo-
rovigatti.github.io/oxDNA/oat/index.html.

##  6 SIMULATION VISUALIZATION

The oxDNA viewer, or oxView, is a robust, browser-based tool for visualizing
and editing oxDNA configurations. Using the Three.js JavaScript library, it
provides a seamless experience for handling large configuration files, with
the capacity to manage up to 1.2 million nucleotides. Its interface is
intuitive, allowing users to simply drag and drop a topology, configuration,
or trajectory files into the browser window. It also supports JSON overlay
files and offers a variety of features, including video options, different
output formats, console commands, and APIs for scene and edit operations.
Additionally, it enables rigid body simulations and 3D printing exports.

Notably, oxView incorporates a "relaxation" simulation, a necessary step prior
to the actual data collection to correct potential physical inconsistencies in
nucleotide placements, which can replace the necessity to run relaxation
simulation from the command line as described above. For a more detailed
understanding of oxView’s capabilities and applications, see [14].

##  7 FUTURE DEVELOPMENT

The simulation code is still in an active development, and as such further
extensions (such as inclusion of more accurate coarse-grained representation
of proteins, nanoparticles, support for DNA-RNA hybrids and more accurate
representation of different salt conditions, sequence-dependent mechanical
properties, capture of hydrodynamic interactions in kinetic simulations) are
either in active development or on the feature-request lists. Similarly,
should more experimental data become available (e.g. thermodynamic
measurements of melting temperatures of large DNA nanostructures, kinetic
measurements of strand displacement reactions and single molecules experiments
on duplex as well as single-stranded DNA or RNA), it will present opportunity
for application of novel fitting methods and verification algorithms to better
parameterize the model.

At the same time, the detailed oxDNA and oxRNA simulations can be used to
parameterize Markov chain models (e.g. [37]) that describe interactions
between DNA and RNA strands in a state space defined by number of base pairs
between strands rather than in terms of x-y-z coordinates in 3D space as
molecular models do. Where limited experimental data are available, the oxDNA
and oxRNA models can be used to augment the training dataset.

Every new release of oxDNA aims to preserve backward compability, even if new
features are added that might make some of the older approaches obsolete. To
get in touch with feature request or ask for help with exisiting code, the
developers can be contacted at https://github.com/lorenzo-
rovigatti/oxDNA/issues.

## ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation
under Grant DMR-2239518. We acknowledge contributions of Thomas Ouldridge,
Flavio Romano, Ben Snodin, Erik Poppleton, and Lorenzo Rovigatti to the oxDNA
model and simulation and analysis tools development.

## References

  * [1] N. C. Seeman.  Nucleic acid junctions and lattices.  Journal of theoretical biology, 99(2):237–247, 1982.
  * [2] N. C. Seeman and H. F. Sleiman.  DNA nanotechnology.  Nature Reviews Materials, 3(1):1–23, 2017.
  * [3] A. V. Pinheiro, D. Han, W. M. Shih, and H. Yan.  Challenges and opportunities for structural DNA nanotechnology.  Nature nanotechnology, 6(12):763–772, 2011.
  * [4] P. W. Rothemund.  Folding DNA to create nanoscale shapes and patterns.  Nature, 440(7082):297–302, 2006.
  * [5] A. Sengar, T. E. Ouldridge, O. Henrich, L. Rovigatti, and P. Šulc.  A primer on the oxDNA model of DNA: when to use it, how to simulate it and how to interpret the results.  Frontiers in Molecular Biosciences, 8:693710, 2021.
  * [6] C. Maffeo and A. Aksimentiev.  MrDNA: a multi-resolution model for predicting the structure and dynamics of DNA systems.  Nucleic acids research, 48(9):5135–5146, 2020.
  * [7] T. E. Ouldridge, A. A. Louis, and J. P. Doye.  Structural, mechanical, and thermodynamic properties of a coarse-grained DNA model.  The Journal of chemical physics, 134(8):02B627, 2011.
  * [8] P. Šulc, F. Romano, T. E. Ouldridge, L. Rovigatti, J. P. Doye, and A. A. Louis.  Sequence-dependent thermodynamics of a coarse-grained DNA model.  The Journal of chemical physics, 137(13):135101, 2012.
  * [9] P. Šulc, F. Romano, T. E. Ouldridge, J. P. Doye, and A. A. Louis.  A nucleotide-level coarse-grained model of RNA.  The Journal of chemical physics, 140(23):06B614_1, 2014.
  * [10] B. E. Snodin, F. Randisi, M. Mosayebi, P. Šulc, J. S. Schreck, F. Romano, T. E. Ouldridge, R. Tsukanov, E. Nir, A. A. Louis, et al.  Introducing improved structural properties and salt dependence into a coarse-grained model of DNA.  The Journal of chemical physics, 142(23):06B613_1, 2015.
  * [11] J. P. Doye, T. E. Ouldridge, A. A. Louis, F. Romano, P. Šulc, C. Matek, B. E. Snodin, L. Rovigatti, J. S. Schreck, R. M. Harrison, et al.  Coarse-graining DNA for simulations of DNA nanotechnology.  Physical Chemistry Chemical Physics, 15(47):20395–20414, 2013.
  * [12] J. SantaLucia Jr.  A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics.  Proceedings of the National Academy of Sciences, 95(4):1460–1465, 1998.
  * [13] D. H. Mathews, J. Sabina, M. Zuker, and D. H. Turner.  Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure.  Journal of molecular biology, 288(5):911–940, 1999.
  * [14] J. Bohlin, M. Matthies, E. Poppleton, J. Procyk, A. Mallya, H. Yan, and P. Šulc.  Design and simulation of DNA, RNA and hybrid protein–nucleic acid nanostructures with oxView.  Nature protocols, 17(8):1762–1788, 2022.
  * [15] J. Procyk, E. Poppleton, and P. Šulc.  Coarse-grained nucleic acid–protein model for hybrid nanotechnology.  Soft Matter, 17(13):3586–3593, 2021.
  * [16] E. Poppleton, M. Matthies, D. Mandal, F. Romano, P. Šulc, and L. Rovigatti.  oxDNA: coarse-grained simulations of nucleic acids made simple.  Journal of Open Source Software, 8(81):4693, 2023.
  * [17] O. Henrich, Y. A. Gutiérrez Fosado, T. Curk, and T. E. Ouldridge.  Coarse-grained simulation of DNA using LAMMPS: An implementation of the oxDNA model and its applications.  The European Physical Journal E, 41:1–16, 2018.
  * [18] J. N. Zadeh, C. D. Steenberg, J. S. Bois, B. R. Wolfe, M. B. Pierce, A. R. Khan, R. M. Dirks, and N. A. Pierce.  NUPACK: Analysis and design of nucleic acid systems.  Journal of computational chemistry, 32(1):170–173, 2011.
  * [19] R. Lorenz, S. H. Bernhart, C. Höner zu Siederdissen, H. Tafer, C. Flamm, P. F. Stadler, and I. L. Hofacker.  ViennaRNA Package 2.0.  Algorithms for molecular biology, 6:1–14, 2011.
  * [20] L. Rovigatti, J. Russo, F. Romano, M. Matthies, L. Kroc, and P. Šulc.  A simple solution to the problem of self-assembling cubic diamond crystals.  Nanoscale, 14(38):14268–14275, 2022.
  * [21] B. E. Snodin, F. Romano, L. Rovigatti, T. E. Ouldridge, A. A. Louis, and J. P. Doye.  Direct simulation of the self-assembly of a small DNA origami.  ACS nano, 10(2):1724–1737, 2016.
  * [22] D. Frenkel and B. Smit.  Understanding molecular simulation: from algorithms to applications, volume 1.  Elsevier, 2001.
  * [23] M. Tuckerman.  Statistical mechanics: theory and molecular simulation.  Oxford university press, 2010.
  * [24] B. Munsky, W. S. Hlavacek, and L. S. Tsimring.  Quantitative biology: theory, computational methods, and models.  MIT Press, 2018.
  * [25] S. M. Douglas, A. H. Marblestone, S. Teerapittayanon, A. Vazquez, G. M. Church, and W. M. Shih.  Rapid prototyping of 3D DNA-origami shapes with caDNAno.  Nucleic acids research, 37(15):5001–5006, 2009.
  * [26] S. Williams, K. Lund, C. Lin, P. Wonka, S. Lindsay, and H. Yan.  Tiamat: a three-dimensional editing tool for complex DNA structures.  In DNA Computing: 14th International Meeting on DNA Computing, DNA 14, Prague, Czech Republic, June 2-9, 2008. Revised Selected Papers 14, pages 90–101. Springer, 2009.
  * [27] C.-M. Huang, A. Kucinic, J. A. Johnson, H.-J. Su, and C. E. Castro.  Integrated computer-aided engineering and design for DNA assemblies.  Nature Materials, 20(9):1264–1271, 2021.
  * [28] D. Doty, B. L. Lee, and T. Stérin.  scadnano: a browser-based, scriptable tool for designing DNA nanostructures.  arXiv preprint arXiv:2005.11841, 2020.
  * [29] A. Suma, E. Poppleton, M. Matthies, P. Šulc, F. Romano, A. A. Louis, J. P. Doye, C. Micheletti, and L. Rovigatti.  Tacoxdna: A user-friendly web server for simulations of complex DNA structures, from single strands to origami.  Journal of computational chemistry, 40(29):2586–2595, 2019.
  * [30] E. Poppleton, R. Romero, A. Mallya, L. Rovigatti, and P. Šulc.  OxDNA. org: a public webserver for coarse-grained simulations of DNA and RNA nanostructures.  Nucleic acids research, 49(W1):W491–W498, 2021.
  * [31] J. Russo, P. Tartaglia, and F. Sciortino.  Reversible gels of patchy particles: role of the valence.  The Journal of chemical physics, 131(1):014504, 2009.
  * [32] S. Whitelam and P. L. Geissler.  Avoiding unphysical kinetic traps in Monte Carlo simulations of strongly attractive particles.  The Journal of chemical physics, 127(15):154101, 2007.
  * [33] N. Srinivas, T. E. Ouldridge, P. Šulc, J. M. Schaeffer, B. Yurke, A. A. Louis, J. P. Doye, and E. Winfree.  On the biophysics and kinetics of toehold-mediated DNA strand displacement.  Nucleic acids research, 41(22):10641–10658, 2013.
  * [34] F. Hong, J. S. Schreck, and P. Šulc.  Understanding DNA interactions in crowded environments with a coarse-grained model.  Nucleic acids research, 48(19):10726–10738, 2020.
  * [35] M. Sample, M. Matthies, and P. Šulc.  Hairygami: Analysis of DNA Nanostructure’s Conformational Change Driven by Functionalizable Overhangs.  arXiv preprint arXiv:2302.09109, 2023.
  * [36] A. Grossfield.  WHAM: The Weighted Histogram Analysis Method.  Version 2.0.11, http://membrane.urmc.rochester.edu/wordpress/?page_id=126.
  * [37] S. Zolaktaf, F. Dannenberg, M. Schmidt, A. Condon, and E. Winfree.  Predicting DNA kinetics with a truncated continuous-time Markov chain method.  Computational Biology and Chemistry, page 107837, 2023.


11institutetext: 1 Technical University of Denmark, Kongens Lyngby, Denmark
11email: {manli<EMAIL_ADDRESS>
2 Slagelse Hospital, Copenhagen, Denmark
3 Region Hovedstaden. Copenhagen, Denmark

# DTU-Net: Learning Topological Similarity for Curvilinear Structure
Segmentation

Manxi Lin 11 Kilian Zepf 11 Anders Nymark Christensen 11 Zahra Bashir 22
Morten Bo Søndergaard Svendsen 33 Martin Tolsgaard 33 Aasa Feragen 🖂 11

###### Abstract

Curvilinear structure segmentation is important in medical imaging,
quantifying structures such as vessels, airways, neurons, or organ boundaries
in 2D slices. Segmentation via pixel-wise classification often fails to
capture the small and low-contrast curvilinear structures. Prior topological
information is typically used to address this problem, often at an expensive
computational cost, and sometimes requiring prior knowledge of the expected
topology.

We present DTU-Net, a data-driven approach to topology-preserving curvilinear
structure segmentation. DTU-Net consists of two sequential, lightweight
U-Nets, dedicated to texture and topology, respectively. While the texture net
makes a coarse prediction using image texture information, the topology net
learns topological information from the coarse prediction by employing a
triplet loss trained to recognize false and missed splits in the structure. We
conduct experiments on a challenging multi-class ultrasound scan segmentation
dataset as well as a well-known retinal imaging dataset. Results show that our
model outperforms existing approaches in both pixel-wise segmentation accuracy
and topological continuity, with no need for prior topological knowledge.

###### Keywords:

Curvilinear segmentation, topology preservation, triplet loss.

##  1 Introduction

Figure 1: Volumetric and curvilinear structures important for standard plane
assessment in ultrasound of fetal head and abdomen, as well as maternal cervix
are shown in (a)-(c), with manual segmentations in (d)-(f). Please zoom for
details.

Curvilinear structures represent thin and long objects in images, and their
segmentation has extensive application in computer vision and medical imaging,
e.g., segmentation of blood vessels [8], neurons [11], or organ boundaries
[7].

Segmenting curvilinear structures is difficult due to their long, thin shapes:
Their connectivity can be completely altered by misclassifying a small number
of pixels. This is particularly challenging in images with low contrast or
resolution, where the structure blends more with the background. As a result,
connectivity errors due to incorrectly splitting curves (we term these "false
splits") and incorrectly connecting them (we term these "missed splits"), are
very common.

In this paper, we consider curvilinear segmentation problems where missed
splits are just as harmful as false splits. This is in contrast with parts of
the existing literature [13, 3], which focuses on avoiding false splits or
reproducing a prescribed topology by connecting components using prior
knowledge. Our focus on missed splits is motivated by image quality assessment
in fetal ultrasound [18]. Here, image quality depends on particular organs
being completely visible in the image, making over-segmentation directly
harmful. Examples of curvilinear structures important for fetal ultrasound
quality assessment are shown in Figure 1.

Figure 2: A globally topologically correct segmentation with local topological
errors.

####  1.0.1 Curvilinear structure segmentation is a local topological problem.

Curvilinear structure segmentation is often described as a topological problem
[16, 3]. Indeed, the topology of a subset of a topological space refers to
those properties of the set that are invariant under continuous deformation.
In curvilinear structure segmentation, the underlying task is to recover a
curve that is visible as an elongated structure in the image. In such
applications, it is often important to maintain the connectivity of the curve
as accurately as possible. This is a topological problem: If $0<a<b<1$, then
(the image of) a generic curve $c\colon[0,1]\to\Omega$ is not topologically
equivalent to (the image of) its restriction to an interval with a gap $(a,b)$
given by $c|([0,a]\cup[b,1])$.

However, we want to emphasize that correct connectivity is a _local_
topological problem. To this end, note that topology is a global phenomenon,
whereas image segmentation is inherently local. While topology, as quantified
by Betti numbers, checks whether you have found the correct number of
components or loops, this does not ensure a topologically correct
segmentation.

To see this, consider the two curvilinear structures illustrated in Fig. 2,
where the top illustrates a true segmentation and the bottom illustrates an
estimated segmentation. These two segmentations are topologically equivalent,
and their Betti numbers would be identical. They would also be rather similar
in pixel-wise losses or measures. However, from the point of view of topology-
preserving segmentation, the error made is considerable: The ground truth loop
has not been detected, and the loop found in the segmentation is not supposed
to be a loop. In other words, the segmentation needs to be topologically
correct in a local sense. As we shall iterate throughout the paper, our
modelling choices, as well as how we validate the topological correctness of
segmented curvilinear structures, are directly linked to this modelling
rationale.

####  1.0.2 To obtain a data-driven curvilinear segmentation algorithm that is
robust in the face of noisy data, we make the following contributions:

  1. 1.

We build a dual-decoder neural network for curvilinear structure segmentation,
generating predictions by softly fusing the two predictions produced by the
two-stage prediction based on image texture and topology, respectively. The
dual predictions give increased stability for low-quality images.

  2. 2.

Utilizing deep contrastive learning, we leverage the latent embedding
similarity to learn topological information of the images in an end-to-end
scheme, which does not rely on any prior information and produces self-
supervised embeddings that respect the curvilinear structure topology. This
data-driven approach to curvilinear structure segmentation is particularly
scalable when faced with complex topological structures.

  3. 3.

We demonstrate our model on a challenging multi-class ultrasound segmentation
task, and the open DRIVE retinal vessel segmentation dataset.

##  2 Related Work

Image segmentation is usually formulated as pixel-wise classification based on
image texture. However, this is not enough to capture complete curvilinear
structures in the image, as small pixel-wise errors can lead to incorrect
connectivity. Early methods detected curvilinear structures by designing
filters that target their geometrical properties, which are hard to detect at
low contrast [4]. With the rise of deep learning, convolutional neural
networks (CNNs), such as the U-net [12], have become the state-of-the-art for
segmentation.

For curvilinear structure segmentation, one branch of research uses persistent
topology to [3, 7, 6] to train networks to explicitly avoid breaking
topological structure. A caveat with such approaches, especially for complex
topologies, is that they require expensive topological computation. Following
a similar, but less expensive route, Shit et al. [16] define a clDice loss
that encourages the segmentations to maintain correct topology based on their
soft skeletons.

Another branch of work [11, 8] designs delicate backbones and schemes to
capture the details of curvilinear features. CS-Net [11] introduces a self-
attention mechanism that deals with multiple curvilinear types in a unified
manner. IterNet [8] concatenates U-Net and multiple mini U-Nets to incorporate
the structural redundancy. Nevertheless, these methods, which simply learn
semantic labels per pixel, often fail to preserve the connectivity of objects.
TR-GAN [2] is a generative adversarial network that improves artery/vein
classification of vessels by improving segmented vessel connectivity. A
topology ranking discriminator along with a contrastive loss encourages the
network to generate segmentation with the expected topology ranking. The
topological features learned by a pre-trained network are only involved in the
computation of cost functions to preserve the image topology. These methods
have in common that they are not designed to handle low-quality images. The
pre-trained features used in [2, 10] are less well suited for such problems,
and we find experimentally that architectures such as [8] struggle with low-
quality images. In this paper, we, instead, explicitly model the types of
errors for which we wish the network to be robust.

A different alternative is given by post hoc approaches, which correct the
broken predictions from an independent segmentation model. Sasaki et al. [13]
propose an auto-encoder that detects and paints the gaps in images
automatically. These models suffer from the clear disadvantage that they can
only repair false splits – they are unable to handle missed splits.

All the previously mentioned models have in common that they aim to preserve
topology and encode image textures in a single encoding stage. We find that in
particular for images with low contrast or resolution, combining these two
tasks at the same stage may make the model hard to tune. Thus, in this paper,
we split the topology-preserving curvilinear segmentation into two simpler
sub-tasks: First learning texture information, and next learning topological
information from the texture-based coarse segmentation.

##  3 Method

###  3.1 Decoupling Topology-preserving Curvilinear Segmentation

Figure 3: Network architecture of DTU-Net. The dotted arrows refer to
operations that are only performed at the training stage. Please zoom in for
details.

We first summarize our proposed Dual-decoder and Topology-aware U-Net (DTU-
Net) architecture shown in Figure 3. Our DTU-Net includes two lightweight mini
U-Nets, with fewer layers than the original U-Net, referred to as the _texture
net_ and the _topology net_. The mini U-Net architecture is the same as that
of the mini U-Nets from IterNet [8]. The texture net follows the learning
scheme of most general segmentation tasks - encoding image texture information
and giving pixel-wise predictions. The coarse prediction from the texture net
is further fed to the topology net for topology preservation. The topology net
is another mini U-Net, where the encoder learns topological features using a
self-supervised triplet loss, and the decoder summarizes the topology in a
binary segmentation map indicating image foreground and -background. The
texture-based and topology-based predictions are further softly fused to
obtain the final segmentation.

Note that in practice, the texture net can be replaced with any pre-trained or
randomly-initialized segmentation network. This means that our model can also
be employed as an easy-to-use plugin to extend to powerful general
segmentation models. Our use of lightweight mini-U-nets is motivated by the
ultrasound image quality assessment task, where we aim for segmentations to
serve as real-time feedback on image quality during scanning, with limited
computational resources.

###  3.2 Topology Net

####  3.2.1 Learning topological embeddings with a triplet loss.

Recent approaches to curvilinear structure segmentation [2, 10] adopt a pre-
trained VGG19 network to extract high-level image topological features, where
contrastive losses penalize segmentations with different topological features
from the ground truth. In these methods, topological information is extracted
by models pre-trained on natural images, which is agnostic to the task. In
addition, during training, only the feature differences are considered instead
of the features themselves. We, instead, integrate contrastive learning into
the latent embedding of a standard segmentation network using self-supervised
learning to learn topological consistency.

Contrastive learning is often used to extract effective instance features for
downstream tasks by producing embeddings that keep similar objects close, and
push dissimilar objects apart. Here, we adopt a triplet loss in the bottleneck
of the topology net, enabling self-supervised learning of image topology
features.

The triplet loss compares an anchor point to a positive and negative example
[15]. The loss function penalizes the difference between the anchor’s distance
to the positive and negative examples, respectively. In our formulation, the
image is passed first through the texture net $\Theta$ to obtain a coarse
segmentation, and further through the encoder $\Psi$ of the topology net,
where the result $\Psi(\Theta(I))$ is used as an _anchor_. The positive
example $\Psi(G)$ is obtained by passing the ground truth segmentation $G$
through the topology net encoder, and the negative example $\Psi(\hat{G})$ is
similarly obtained from a segmentation $\hat{G}$, which has been corrupted to
alter its topology. As detailed in the next section, we generate positive and
negative examples from segmentation ground truth online during training and do
not need additional annotations.

Sequentially, the triplet loss is defined as:

$L_{tri}(\Psi(\Theta(I)),\Psi(G),\Psi(\hat{G}))=max(d^{+}-d^{-}+\tau,0)$ (1)

where $d^{+}$ is the averaged Euclidean distance of $\Psi(\Theta(I))$ and
$\Psi(G)$, $d^{-}$ is the distance of $\Psi(\Theta(I))$ and $\Psi(\hat{G})$
and $\tau$ is an empirically-set margin. By minimizing the triplet loss during
training, the network forces the encoded features closer to those of the
ground truth segmentation mask whose topology is correct, while pushing it
further away from the embedding of the corrupted segmentation masks.

####  3.2.2 Self-supervised training of the triplet loss.

We generate corrupted masks from the ground truth mask by randomly breaking
the image topology. Synthetic false splits in the corrupted masks are
generated as follows: We assume that false splits can happen anywhere on the
curves when the CNNs fail to recognize the objects. Therefore, segmentation
labels are split into small patches, randomly removing $\lambda$ of the
patches containing foreground information, where $\lambda\in[0,1]$ is a hyper
parameter controlling the degree of corruption.

Next, we also generate corrupted masks with synthetic missed splits: We deem
the missed splits tend to appear in regions with similar texture features, put
simply in raw images, the pixel values. We compute the average $m$ and
standard deviation $d$ of the foreground pixel values in the image (foreground
pixels are given in the ground truth). To generate synthetic missed splits, a
fixed percentage $\lambda$ of the pixels with value within $[m-d,m+d]$ are
selected randomly and assigned the label of the most similar foreground pixel.

In our implementation, adding missed splits and false splits is done
sequentially. The introduction of missed splits brings not only incorrect
connections between curves but also noise, which harms the topology as well.

During inference, the encoder only forwards $\Psi(\Theta(I))$. The triplet
loss in Eq. 1 quantifies the pixel-wise feature difference. We also merge the
global topological feature into the embedding by a SENet [5], which takes the
average of the learned feature map and encodes them into attentive scores via
a multi-layer perceptron. The scores are then multiplied by the original
feature map.

####  3.2.3 Fine-to-coarse topology summarization.

In the DTU-Net, the texture net learns low-level information directly from raw
images, which is effective in semantic segmentation tasks. Our topology net,
however, abstracts the information into high-level topological features, which
are less applicable for texture recognition tasks due to the loss of texture
details. We expect the topology net to summarize the information from object
connectivity, and contribute to the final segmentation. Consequently, we ask
the topology net a binarized question: Based on the learned topology, should
the structure be connected or not? For multi-class segmentation tasks such as
in our ULTRASOUND dataset, these binary predictions are sequentially employed
to refine the multi-class semantic segmentation from the texture net.

###  3.3 Segmentation Refinement and Joint Feature Learning

Our DTU-Net consists of a texture network and a topology network. We denote
the texture network (a mini U-Net), the encoder, and the decoder of the
topology network (another mini U-Net) as $\Theta(\cdot)$, $\Psi(\cdot)$, and
$\Omega(\cdot)$ respectively. Each part of the network can be trained
separately, as with the post hoc methods. We train the network in a joint way,
where the gradient backpropagation from the topology net can also update the
parameters in the texture net. Denoting input images $I$ and segmentation
ground truth masks $G$, DTU-Net is trained by a unified loss

$L_{DTU}(I,G)=L_{tex}(\Theta(I),G)+L_{BCE}(\Omega(\Psi(\Theta(I))),\bar{G})+L_{tri}$
(2)

where $L_{tex}$ is the pixel-wise segmentation loss for the texture net,
$L_{BCE}$ is the binary cross entropy loss tuning the prediction from the
topology net, $L_{tri}$ is the triplet loss in Eq. 1, $\bar{G}$ and $\hat{G}$
refers to the ground truth binarized and corrupted by us. The loss function
also indicates that our DTU-Net does not need additional annotations and
learns image topology in a self-supervised way.

As the topology net and the texture net view the segmentation problem from
different points of view, we fuse their predictions into a final prediction –
this avoids sensitivity to the topology net sometimes being overconfident.
Inspired by ensemble learning [14], we fuse the two predictions using weighted
summation, where the hyper parameter $\omega$ indicates our confidence in the
texture-based prediction. For each pixel, given $p_{top}$ the pixel-wise
probability of the topological prediction, and $p_{tex}^{(i)}$ the pixel-wise
probability of the $i$-th category from the textural prediction, the two
predictions are simply fused by:

$p^{(i)}_{final}=\begin{cases}(1-\omega)(1-p_{top})+\omega
p^{(0)}_{tex},i=0;\\\
\frac{(1-\omega)p_{top}+\omega\sum_{i=1}^{c}p^{(i)}_{tex}}{\sum_{i=1}^{c}p^{(i)}_{tex}}p^{(i)}_{tex},i=1,2,...,c\end{cases}$
(3)

where $c$ is the total number of categories in the semantic segmentation task,
$p^{(i)}_{final}$ is the pixel-wise final prediction, the $0$-th category
means the background.

##  4 Experiments

####  4.0.1 Baseline models.

We compare our model against a mini U-net trained with the clDice loss
function, as well as with the persistent homology-based TopoLoss [7]; the
data-driven IterNet [8] and the inpainting network of [13]. The input of the
inpainting network is the prediction from a mini U-Net trained by a focal
loss.

The clDice loss function preserves topological information by learning object
skeletons, which, however, does not guarantee accurate segmentation. In the
original paper [16], this is handled by mixing the clDice metric with another
pixel-wise classification loss, more precisely the soft dice loss with hyper-
parameter $\beta$:

$L_{clDice}=\beta(1-clDice)+(1-\beta)(1-softDice)$ (4)

In our experiments, and in particular on the more challenging ULTRASOUND
dataset, we found the clDice loss function challenging to tune. The
experimental setup from [16] did not converge with of the suggested values of
$\beta=0$ or $0.5$ (see Table3(a)). To improve convergence, we replaced
softDice with focal loss [9]:

$L_{clDice}=\alpha(1-clDice)+(1-\alpha)L_{focal}.$ (5)

Even here, however, the performance was unstable with the choice of $\alpha$,
as can be seen from Table 3(a), and we choose to compare to the optimal
observed $\alpha=0.2$ in Table 1. For the DRIVE dataset, we choose
$\alpha=0.5$ following [16].

Table 1: Results on the ULTRASOUND dataset. Note that cIoU and cmIoU refer to IoU and mean class IoU, respectively, for curvilinear structures, whereas vIoU and vmIoU refers to volumetric structures. Model | Loss | Frechet | Betti error | vIoU | vmIoU | cIoU | cmIoU
---|---|---|---|---|---|---|---
mini U-Net | clDice | 0.7801 | 0.2635 | 97.95 | 29.84 | 31.41 | 20.41
mini U-Net | TopoLoss | 0.9653 | 0.2606 | 98.07 | 26.29 | 25.89 | 18.56
IterNet | focal | 0.4329 | 1.0963 | 98.20 | 17.67 | 5.76 | 3.70
inpainting | MSE | 0.8189 | 1.1179 | 93.88 | 23.93 | 40.15 | 24.37
DTU-Net | Eq. 2 | 0.5818 | 0.0974 | 99.51 | 44.33 | 72.30 | 52.92
Table 2: Segmentation results on the DRIVE dataset. Model | Loss | Frechet | Betti error | IoU
---|---|---|---|---
mini U-Net | clDice | 3.6850 | 1.0890 | 67.95
mini U-Net | TopoLoss | 3.6070 | 0.9232 | 70.52
IterNet | focal | 3.5335 | 1.0039 | 70.35
inpainting | MSE | 3.9001 | 1.3824 | 65.54
DTU-Net | Eq. 2 | 2.9316 | 0.8597 | 73.86

####  4.0.2 Tasks.

We validate our algorithm on two different segmentation tasks.

Table 3: (a) clDice model selection on ULTRASOUND. The first two rows show the
models suggested in the original paper [16]; these did not converge. The next
rows show results for various parameters using the focal loss in place of the
soft dice. (b) Training times on ULTRASOUND. TopoLoss and clDice were used
with a mini U-Net; the inpainting included two stages.

(a) $\alpha$ | $\beta$ | Frechet | Betti error | vIoU | vmIoU | cIoU | cmIoU
---|---|---|---|---|---|---|---
0 | 0.5 | - | 9.8019 | 2.14 | 0.33 | 1.53 | 1.21
0 | 1 | - | 13.93 | 3.72 | 0.62 | 2.91 | 1.80
0 | 0 | 0.8418 | 0.4572 | 97.70 | 29.42 | 32.99 | 23.14
0.1 | 0 | 0.9500 | 0.3800 | 97.74 | 28.13 | 26.34 | 17.72
0.2 | 0 | 0.7801 | 0.2635 | 97.95 | 29.84 | 31.41 | 20.41
0.3 | 0 | 1.0798 | 0.4128 | 97.70 | 28.75 | 29.79 | 21.07
0.4 | 0 | 1.0727 | 0.3873 | 97.44 | 29.30 | 25.56 | 18.09
0.5 | 0 | 0.8969 | 0.3264 | 97.06 | 26.69 | 22.76 | 15.00

(b) Model | Time/hours
---|---
TopoLoss | 43.74
IterNet | 20.02
clDice | 6.91
inpainting | 7.47
DTU-Net | 17.2

The ULTRASOUND dataset poses a challenging fetal ultrasound multiclass
segmentation task including both volumetric and curvilinear structures. These
images, collected from the Danish national fetal ultrasound screening
database, are obtained a typical 3rd-trimester growth scans, where four
standard ultrasound planes (head, abdomen, femur, cervix) are acquired to
assess preterm birth and fetal weight. Whether the images are "standard"
depends on the visibility and pose of certain anatomical organs, several of
which appear as curvilinear structures in the image. It is important to
accurately detect whether the entire structure is present, as well as
measuring its dimensions. Our motivating downstream application requires
segmenting multiple different organs in real time while scanning; we have
trained our models to segment both the volumetric and curvilinear structures.
There are 2088 images of resolution $960\times 720$, divided into $80\%$ for
training, $10\%$ for validation, and $10\%$ for testing. The dataset includes
6 curvilinear (cervix canal outline, cervix outer boundary, cervix inner
boundary, outer skin boundary, outer bone boundary, thalamus) and 7 volumetric
segmentation categories (bladder. femur bone, stomach bubble, umbilical vein,
kidney, fossa posterior, cavum septi pellucidi). There are 468 femur, 531
abdomen, 686 head, and 403 cervix images. Images are reshaped to $224\times
288$.

The DRIVE retinal segmentation dataset [17] induces a binary segmentation
task. We used $16$ images for training and $4$ for testing. During training,
the images are randomly cropped into $256\times 256$ patches, while in testing
patches are obtained with a $128\times 128$ window sliding with a stride of
$64$.

####  4.0.3 Experimental settings.

The texture net is trained with a combined dice and focal loss as $L_{tex}$ in
Eq. 2. Triplet loss examples are generated incrementally by reducing $\lambda$
from 50% to 10% during training, with $\omega=0.5$ and $\tau=0.1$.

####  4.0.4 Evaluation metrics.

To measure geometric affinity between curves, we employ the Frechet distance
[1], which measures the extent to which the two curves remain close to each
other throughout their course. We find this more suitable than e.g. accuracy,
which is sensitive to small shifts or differences in curve thickness.

To measure the local topological correctness of the segmented structures, we
apply a sliding window, in which we compare the number of components in the
prediction and the ground truth annotation, respectively. This corresponds to
a local Betti number computation. For each image, we compute a topological
error (Betti error) which is the average difference in Betti numbers between
ground truth and prediction over all patches seen, similar to what is done in
[7].

For completeness, we also include the voxel-wise intersection over union
(IoU), and mean class IoU for the multi-class segmentation task. For the
ULTRASOUND dataset we report these both for the curvilinear and volumetric
structure classes.

Table 4: Ablation study on different components of DTU-Net. Model | Frechet | Betti error | vIoU | vmIoU | cIoU | cmIoU
---|---|---|---|---|---|---
w/o triplet loss | 0.6605 | 0.1056 | 98.95 | 31.56 | 66.45 | 41.21
$\omega=0$ | 0.5863 | 0.1023 | 98.90 | 43.11 | 70.43 | 52.33
$\omega=1$ | 0.8601 | 0.1683 | 98.12 | 14.15 | 30.88 | 22.20
DTU-Net | 0.5818 | 0.0974 | 99.51 | 44.33 | 72.30 | 52.92
Table 5: Ablation on Model size Model | Size | Frechet | Betti error | vIoU | vmIoU | cIoU | cmIoU
---|---|---|---|---|---|---|---
full U-Net | 131.8M | 0.5849 | 0.1383 | 98.53 | 31.09 | 32.58 | 23.22
DTU-Net | 65.5M | 0.5818 | 0.0974 | 99.51 | 44.33 | 72.30 | 52.92

####  4.0.5 Results.

Results for the ULTRASOUND dataset are found in Table 1, and for the DRIVE
dataset in Table 2, with example segmentations in Figure 4. The performance of
the model used as the input for the post hoc inpainting on the ULTRASOUND
dataset can be found in Table 3(a) with $\alpha$ and $\beta=0$.

####  4.0.6 Running time.

Table 3(b) compares training times on ULTRASOUND. The models were trained on a
Quadro RTX 6000 GPU with identical training settings. Note how DTU-Net has a
considerably lower training time than TopoLoss, which is the closest
competitor on the DRIVE dataset. While clDice and inpainting have the lowest
runtimes, these do not compete in performance (Tables 1 and 2).

####  4.0.7 Ablation study.

We conducted an ablation study over different components of the DTU-Net on the
ULTRASOUND dataset, see Table 4. We see that removing the triplet loss from Eq
2 leads to increased Betti error. We also evaluate the effectiveness of
feature fusion. Specifically, $\omega$ in Eq 3 is set to be 0 and 1
respectively, meaning that the prediction of foreground and background is
dominated by the topology net or the texture net. To prove that our model
performance is not merely due to the increase of learnable parameters, we also
include a comparison with a larger U-Net trained with the clDice loss
($\alpha=0.2$), which has roughly two times the model size of ours (Table 5).

##  5 Discussion and Conclusion

We have designed a framework for curvilinear structure segmentation which is
trained, in a data driven way, to avoid both false and missed splits. Our
method is built around a dual-decoder network, which first makes a rough
segmentation based on image texture, and next refines this optimizing for
predicting correct topological structures. Here "correctness" also emphasises
whether the structure is visible in the image. As the second network takes the
predicted foreground probabilities as input, it learns to enhance those
responses that form curvilinear structures. This is obtained using self-
supervised contrastive learning, to teach the network the difference between a
weak curvilinear response, and no response. On the highly challenging
ULTRASOUND dataset, we see that our method outperforms all other methods. This
is natural as most methods are developed using easier segmentation tasks such
as e.g. vessel segmentation, and they are not optimized to handle the
challenging nature of ultrasound data.

Figure 4: Segmentation of examples from ULTRASOUND and DRIVE. Blue arrows
point at missed splits while red arrows point at false splits in the image.

However, our method also outperforms baselines on the well-known DRIVE
dataset. Here, we hypothesize that the data driven approach may have an easier
time learning to handle complex topologies. It also has an advantage over
persistent homology based methods, as it explicitly seeks to preserve local
topology, whereas persistent homology based methods achieve locality by
preserving Betti numbers within patches. This means that the locality scale
becomes a parameter, which is additionally constrained by computational
demands. This gives the DTU-net an advantage both in terms of choice of scale
(it is inherently robust to scale) and in terms of computational time, as also
shown by our experiments.

One important difference between methods is the number of parameters used. We
have designed most methods around the fixed size mini U-net. Since our own
method has the additional topology network, it de facto has more parameters
than the mini U-net trained with clDice. To asses whether this difference in
parameters was unfair, we also compared our method to a deeper U-net (with
more parameters than our full architecture), trained with clDice. However, the
results (see Table 5) show that we still outperform clDice on the ULTRASOUND
dataset even when it is given the advantage of more parameters. The IterNet
architecture also has more parameters than ours, but is not competitive.

Limitations. Our method is naturally limited by our downstream need for real-
time inference, which is why we rely on mini U-Nets rather than more complex
models. Allowing larger models might change the outcomes for competing
methods, although this was not the case for clDice on the ULTRASOUND dataset.
At the same time, our architecture’s two mini U-Nets makes our method more
expensive, also at inference time, than single mini U-Nets trained with
topological loss functions. Finally, while the data-driven contrastive
learning is clearly helpful in learning topological structure, it does not
come with theoretical guarantees.

Summary. We have proposed a data-driven framework for segmenting curvilinear
structures, and shown that it outperforms the baselines on DRIVE, and by far
on the challenging ULTRASOUND dataset. We find it interesting that all
baseline models struggle on the ultrasound data, and see this as a cue for the
community to include more challenging segmentation tasks when developing
methods for curvilinear structure segmentation.

## References

  * [1] Alt, H., Godau, M.: Computing the Fréchet distance between two polygonal curves. International J. Computational Geometry & Applications 5(01n02), 75–91 (1995)
  * [2] Chen, W., Yu, S., Wu, J., Ma, K., Bian, C., Chu, C., Shen, L., Zheng, Y.: TR-GAN: topology ranking GAN with triplet loss for retinal artery/vein classification. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 616–625. Springer (2020)
  * [3] Clough, J., Byrne, N., Oksuz, I., Zimmer, V.A., Schnabel, J.A., King, A.: A topological loss function for deep-learning based image segmentation using persistent homology. IEEE Transactions on Pattern Analysis and Machine Intelligence (2020)
  * [4] Frangi, A.F., Niessen, W.J., Vincken, K.L., Viergever, M.A.: Multiscale vessel enhancement filtering. In: International conference on medical image computing and computer-assisted intervention. pp. 130–137. Springer (1998)
  * [5] Hu, J., Shen, L., Sun, G.: Squeeze-and-excitation networks. In: IEEE CVPR (2018)
  * [6] Hu, X., Wang, Y., Li, F., Samaras, D., Chen, C.: Topology-aware segmentation using discrete morse theory. In: ICLR (2021)
  * [7] Hu, X., Li, F., Samaras, D., Chen, C.: Topology-preserving deep image segmentation. Advances in Neural Information Processing Systems 32 (2019)
  * [8] Li, L., Verma, M., Nakashima, Y., Nagahara, H., Kawasaki, R.: Iternet: Retinal image segmentation utilizing structural redundancy in vessel networks. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 3656–3665 (2020)
  * [9] Lin, T.Y., Goyal, P., Girshick, R., He, K., Dollár, P.: Focal loss for dense object detection. In: IEEE ICCV (2017)
  * [10] Mosinska, A., Marquez-Neila, P., Koziński, M., Fua, P.: Beyond the pixel-wise loss for topology-aware delineation. In: IEEE CVPR (2018)
  * [11] Mou, L., Zhao, Y., Fu, H., Liu, Y., Cheng, J., Zheng, Y., Su, P., Yang, J., Chen, L., Frangi, A.F., et al.: Cs2-net: Deep learning segmentation of curvilinear structures in medical imaging. Medical image analysis 67, 101874 (2021)
  * [12] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [13] Sasaki, K., Iizuka, S., Simo-Serra, E., Ishikawa, H.: Joint gap detection and inpainting of line drawings. In: IEEE CVPR (2017)
  * [14] Schapire, R.E.: A brief introduction to boosting. In: IJCAI. vol. 99, pp. 1401–1406. Citeseer (1999)
  * [15] Schroff, F., Kalenichenko, D., Philbin, J.: Facenet: A unified embedding for face recognition and clustering. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 815–823 (2015)
  * [16] Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P., Bauer, U., Menze, B.H.: clDice-a novel topology-preserving loss function for tubular structure segmentation. In: IEEE CVPR (2021)
  * [17] Staal, J., Abràmoff, M.D., Niemeijer, M., Viergever, M.A., Van Ginneken, B.: Ridge-based vessel segmentation in color images of the retina. IEEE Transactions on Medical Imaging 23(4), 501–509 (2004)
  * [18] Wu, L., Cheng, J.Z., Li, S., Lei, B., Wang, T., Ni, D.: Fuiqa: fetal ultrasound image quality assessment with deep convolutional networks. IEEE transactions on cybernetics 47(5), 1336–1349 (2017)