11institutetext: Carnegie Mellon University 22institutetext:
AerospaceComputing 33institutetext: NASA Ames 44institutetext: IPSSE, TU
Clausthal

# Verifying the Safety of a Flight-Critical System

Guillaume Brat 11 David Bushnell 22 Misty Davies 33
Dimitra Giannakopoulou 33 Falk Howar F. Howar did this work while at Carnegie
Mellon University.44 Temesghen Kahsai 11

###### Abstract

This paper describes our work on demonstrating verification technologies on a
flight-critical system of realistic functionality, size, and complexity. Our
work targeted a commercial aircraft control system named Transport Class Model
(TCM), and involved several stages: formalizing and disambiguating
requirements in collaboration with domain experts; processing models for their
use by formal verification tools; applying compositional techniques at the
architectural and component level to scale verification. Performed in the
context of a major NASA milestone, this study of formal verification in
practice is one of the most challenging that our group has performed, and it
took several person months to complete it. This paper describes the
methodology that we followed and the lessons that we learned.

##  1 Introduction

This paper demonstrates the use of verification approaches on a safety-
critical system of realistic functionality, size, and complexity. The work
addresses a major milestone of the NASA Aviation Safety program and was
performed over several months by a team involving four verification experts,
one senior software engineer, and an aerospace engineer.

The target of our study is a Simulink model of a twin-engine aircraft
simulation named Transport Class Model (TCM). The TCM was selected for a
number of reasons. First, it is unclassified and can therefore be shared
outside of NASA. This is important to us because we would like the community
to benefit from our experience and to be able to use this as a common
benchmark on which additional verification technologies can be applied.
Second, the system was developed independently by a different NASA center and
therefore we had no prior knowledge of its potential errors, or its design.
The setting was therefore similar to one in which a safety-critical system is
handed to verification experts for analysis and certification, where the
experts were not involved in the system design.

As the TCM does not come with requirements, we used several sources such as
pilot training manuals and the Federal Aviation Regulations for commercial
aircraft, to collect relevant requirements for our study. A significant amount
of our work involved formalizing and disambiguating requirements in
collaboration with domain experts. The resulting requirements constitute
verification properties, which we encoded as synchronous observers in
Simulink.

The Simulink models that describe the system had to be processed in order to
be usable by the verification tools utilized in this study. We used SMT-based
model checking, and in particular the PKind [22] to verify the properties. The
Simulink models including the synchronous observers were automatically
translated into the synchronous dataflow language Lustre [6], to be processed
by PKind.

A major goal of this work was to experiment with compositional techniques to
enable the scalable use of formal methods for systems of realistic size.
Compositional verification constructs a verification argument for a complex
system by composing simpler verification results at the level of the system
components. There are several well-known advantages to taking such an
approach. Scalability of verification is a major driver. Through the
decomposition of system-level requirements into component-level ones at design
time, it is easier to assign clear responsibilities to the developer of each
component. Finally, as components of a system change or evolve, compositional
verification enables the reuse of verification results of unchanged
components.

The work described in this paper is one of the most challenging verification
exercises that our team has performed. As such, it forced us to define a high-
level methodology for the verification of flight-critical systems, and has
enabled us to comment on advantages and limitations of verification techniques
and tools in handling such systems. We found that close collaboration between
verification and domain experts is required to formalize requirements and, in
particular, assumptions about the physical system, without which verification
would fail. Compositional verification is key in constructing scalable and
meaningful proofs for complex systems in the aeronautics domain.

The remainder of the paper is organized as follows. Section 2 describes the
TCM, while section 3 discusses the process through which we obtained
requirements and disambiguated and formalized them for verification. The
verification effort itself is presented in section 4. The experience and
lessons learned from this substantial effort are discussed in section 5, with
section 6 placing this effort in the context of related work and section 7
concluding the paper and discussing future work.

##  2 The Transport Class Model

Our target system is derived from NASA Langley’s Transport Class Model (TCM)
[20], a simulator of a mid-size (approximately 250,000 lb.), twin-engine,
commercial transport-class aircraft. The TCM is not intended as a high-
fidelity simulation of any particular transport aircraft. Rather, it is meant
to be representative of the types of non-linear behaviors of this class of
aircraft.

The TCM includes models for the avionics (with transport delay), actuators,
engines, landing gear, nonlinear aerodynamics, sensors (including noise),
aircraft parameters, equations of motion, and gravity. It is primarily
implemented in Simulink, consisting of approximately 5700 Simulink blocks. The
system also includes several thousands of lines of C/C++ code in libraries,
primarily used for the engines and the nonlinear aerodynamics models. Our work
studies the guidance and controls models and their properties, within the
overall context of the TCM system, and these are implemented entirely in
Simulink.

Figure 1 depicts the top-level controls Simulink (data-flow) diagram. An
aircraft can be controlled either manually (using a control stick and pedals),
or through the mode control panel (MCP) of the flight computer (autopilot). In
the diagram, pink-shaded boxes with dark outlines (to the left) highlight the
inputs from the pilot and the MCP. Red arrows identify the various controls
subsystems (each subsystem may itself be a complex collection of subsystems),
explained below. The TCM contains inner loop proportional-integral (PI)
controllers for all three angular axes of motion (roll, pitch, and yaw). These
inner loop controllers function regardless of whether the pilot flies manually
or the autopilot issues commands. Additionally, the TCM’s autopilot can
control altitude (either by flying to a directed altitude, or holding a
current altitude), can reach and maintain a desired flight path angle (FPA),
can reach and maintain a desired heading, and can control the airplane’s
speed. Finally, the blue-shaded box with the dark outline (to the left) in the
diagram shows the collected outputs—the commands to the actuators.

Figure 1: The Simulink guidance and controls system for the TCM.

##  3 Requirements: Elicitation and Formalization

Written requirements are not available for the autopilot and controls software
of the TCM because it was not intended for embedded production-level code. The
original implementers of the software informed us that their release process
was based on a side-by-side comparison of behaviors between the simulator and
an experimental aircraft. However, for the purpose of our study, we need
safety properties representative of those used in the certification of civil
aviation transport vehicles. For this reason, NASA’s Armstrong Flight Research
Center chose relevant requirements from the Federal Aviation Regulations
(FARs) governing commercial aviation transport vehicles (Part 25) [10], such
as the following:

> _FAR-25.672b_ : The design of the stability augmentation system or of any
> other automatic or power operated system must permit initial counteraction
> of failures of the type specified in section 25.671(c) without requiring
> exceptional pilot skill or strength, by either the deactivation of the
> system, or a failed portion thereof, or by overriding the failure by
> movement of the flight controls in the normal sense.

FAR requirements such as the one above refer to several high-level systems,
but in our study we focus on the Guidance, Navigation, and Control (GNC)
system of Figure 1. The GNC is divided into three types of components: (i)
_mode logic_ capturing the modes in which the autopilot can operate, how these
modes are enabled and disabled, and which control systems are enabled by
particular modes, (ii) components targeting _controllability_ of the airplane
by ensuring that the actuators can, at any instant, respond appropriately to a
command, and (iii) components ensuring the _stability_ and _maneuverability_
of the aircraft, making the aircraft robust to state disturbances and easily
controlled by the pilot.

The TCM does not contain detailed Simulink models for Navigation. Moreover,
classical mathematical techniques like Lyapunov theory [23] provide well-
understood ways of checking stability and maneuverability properties. The
formal verification tools that our study targets are better suited for
checking mode logic and simple controllability properties of type (i) and (ii)
above. As a consequence, we focus on GNC-level requirements, and sub-
requirements (also called “child"-requirements) that refer exclusively to
Guidance. A child-requirement of the above FAR requirement, is, for example,
the following:

> _GNC-150_ : The Guidance Navigation and Control Function shall enable the
> pilot to transition the vehicle from one flight condition to another (i.e.
> climb to level flight) under all operating conditions including failure of a
> single engine.

In order to elicit such higher-level requirements into properties that can be
checked on the components of our case-study, we examined pilot training
materials for the Boeing 737 (B737) Automatic Flight Systems [27]. B737 is
within the class of vehicles that the TCM simulator targets. Within these
training materials were behavioral specifications for the B737 Altitude
Acquire, Altitude Hold, and Level Change Modes, equivalent in functionality to
the TCM’s Altitude Controller. We also used the specifications for the B737
Heading Select Mode as the desired properties for the TCM Heading Controller,
the B737 Glide Slope Capture mode for the TCM Flight Path Angle Controller,
and the B737 Autothrottle for the TCM’s Autothrottle. We started with 88 GNC
requirements, which we tried, with the help of these B737 documents, to map
into safety properties of the Guidance system. This effort resulted into 20
properties, illustrated in Table 1. These 20 properties are sub-requirements
of GNC-level requirements, such as the GNC-150 requirement shown above. In the
rest of this paper, we will use property G-120 as a running example:

> _G-120_ : The guidance shall be capable of climbing at a defined rate, to be
> limited by minimum and maximum engine performance and airspeeds.

Table 1: A summary of verified properties on the TCM. # | Property | Assumptions | Original Requirement
---|---|---|---
1 | G-250 | G-260 | The heading control mode, when selected, sends roll
|  |  | commands to turn to and maintain the commanded heading.
2 | G-110 | G-220,G-260 | The guidance system shall be capable of steering to and
|  |  | following a specified heading.
3 | G-120 | G-180,A1,A2, | The guidance shall be capable of climbing at a defined rate, to be
|  | FPA1 | limited by minimum and maximum engine performance and airspeeds.
4 | G-130 | G-180,A1,A2 | The guidance shall be capable of descending at a defined rate, to be
|  |  | limited by minimum and maximum engine performance airspeeds.
5 | G-140 | G-120,G-200 | The guidance shall be capable of climbing at a specified rate
|  |  | to a specified altitude, to be limited by maximum engine
|  |  | performance for a set airspeed
6 | G-150 | G-180,A1,G-120, | The guidance shall be capable of descending at a specified rate
|  | A2,G-200 | to a specified altitude, to be limited by maximum engine
|  |  | performance for a set airspeed
7 | G-170 (Mode) | – | The altitude control shall engage when the altitude control mode
|  |  | is selected and when the FPA control mode is not selected, and when
|  |  | there is no manual pitch or manual roll command from the stick.
8 | G-180 (Mode) | – | The FPA control shall engage when the FPA mode
|  |  | is selected, and when there is no manual pitch or manual roll
|  |  | command from the stick.
9 | G-100 | – | The Guidance system shall be capable of maintaining a steady speed
|  |  | in the normal flight envelope.
10 | G-200 | – | If the altitude control is engaged, once the plane is within 250 ft of
|  |  | the commanded altitude, the plane will remain within 250 ft
|  |  | of the commanded altitude.
11 | G-210 (Mode) | – | If the FPA control and the altitude control are both selected, the FPA
|  |  | control will disengage and the altitude control will engage once the
|  |  | lane is within 200 ft of the commanded altitude.
12 | G-220 (Mode) | – | The heading control shall engage when the heading control mode
|  |  | is selected, and when there is no manual pitch or manual roll
|  |  | command from the stick.
13 | G-230 | – | If the altitude control is engaged with no active speed control,
|  |  | the speed control shall engage and the speed command shall synchronize
|  |  | to the current speed, which shall become the new altitude’s target speed.
14 | G-240 | – | The bank angle limit is established by the Bank Angle Limit Selector.
15 | G-260 (Mode) | – | When the heading control mode is engaged, roll commands
|  |  | are given to turn in the nearest direction to the selected heading.
16 | G-270 (Mode) | – | Manually positioning the thrust levers does not cause
|  |  | autothrottle disengagement.
17 | G-290 | – | The autothrottle will be limited by the max and the min throttle.
18 | G-160 | – | The guidance function shall be able to automatically deploy spoilers
|  |  | to limit speed in a descent, or when a significant reduction in
|  |  | airspeed is requested by the pilot, deactivating at low speed.
19 | G-280 | – | The FCCs shall issue a warning when the commanded altitude
|  |  | disagrees with the stored commanded altitude stored in the FCCs.
20 | G-190 | – | If any control surface actuator loses hydraulic pressure,
|  |  | the autopilot shall disengage.

###  3.1 Formalization

For verification, the properties shown in Table 1 must be disambiguated and
formulated in terms of the signals (inputs and outputs) of the TCM model.
Moreover, they have to be written in a formal language which, in our case, is
Simulink (see Section 4). The requirements formalization process was performed
over several iterations, and involved discussions between domain and
verification experts. It included several steps, presented in this section.

1\. Develop a shared understanding of the requirement. Natural language often
allows for slightly different interpretations. In some cases, properties even
interfered or contradicted one another. For example, what does “… to be
limited by minimum and maximum engine performance …” mean in the context of
G-120? It could mean that a. the guidance shall be capable to climb when
commanded to do so, but the climb rate should be limited by engine performance
and airspeeds, or that b. the guidance is only required to be able to climb
when the defined rate is within the minimum and maximum engine performance.
For all such cases, we consulted with the domain expert on the team to develop
a common understanding and/or we refined properties to be more precise (e.g.,
interpretation a. was selected for property G-120).

2\. Decompose the property into a requirement on the control system and
assumptions about the physics of the airplane. Since we reason formally only
about controllability, we decomposed the properties into a part that can be
proven formally and into corresponding assumptions about stability and
physics. Such assumptions were not formally verified, but were confirmed by
domain experts. For example, while G-120 specifies that the aircraft shall
climb at a defined rate, the control system has got just sensors and
actuators. Since we do not mathematically specify the physics behavior, we do
not have a formal definition of what it means to climb at a defined rate in
terms of the sensor values—instead we define climbing solely in terms of the
actuator commands. Since we focus on instantaneous controllability, we require
that the control system outputs a value to an actuator that moves the aircraft
into the right direction (e.g., if the current climb rate is smaller than the
commanded climb rate, the control system should issue a command to the
ailerons that would pitch the aircraft upwards). We then assume that the inner
loop controllers and physics will result in an increased climb rate.

3\. Identify affected components and signals of the control model. At this
stage properties are still formulated in natural language and state vague
things like “The guidance shall be capable …”. Such expressions cannot be
mapped to the TCM directly. The model describes sensors and signals. We
therefore have to express properties in terms of the signals available in our
model. In our study, we motsly relied on domain experts to help us with this
step. In the case of G-120, we had to refine “being capable” as the concrete
situations in which the control system is expected to act. This could be
expressed as “If in FPA-control mode, and if there is no manual aileron or
pitch command from the pilot …”, which can be mapped to signals in the model.

4\. Decompose the requirement on the control system into sub-requirements on
single components. In some cases, the formalized requirements specify behavior
of the complete control system in terms of its global inputs and outputs.
However, proving the requirements marked gray (7-17) in Table 1 required
information about internal signals between lower level components. We
therefore decomposed these requirements into sub-requirements over internal
signals between components. These sub-requirements were not merely slices of
the global property but actual assume/guarantee pairs that we derived
manually.

Property G-120, for example, was decomposed as illustrated in Fig. 3 (the
figure displays Lustre code, as translated from Simulink by our compiler). The
decomposition expresses the fact that the FPA control module, when engaged, is
in charge of maintaining an FPA (FPA1). The mode logic ensures that the
guidance system cannot be in Altitude Control mode and FPA control mode at the
same time (i.e., these modes are mutually exclusive). Based on this fact, the
remaining properties express that: – in FPA control mode, the FPA control
module is engaged (G-180); – if not in Altitude control mode, then the
Altitude control module is not engaged (A1), and when not engaged, the
Altitude control module will not send commands (A2). As a result, when in FPA
control mode, the FPA will be the only mode engaged, and it will issue
commands to maintain an FPA, which means that the guidance system is able to
climb at a defined rate.

Property formalization helped us identify missing functionalities in the TCM
model and also led us to refine many properties.

  * •

Three requirements specify behavior for components not modeled by the TCM
(e.g., spoilers). These requirements could not be formalized.

  * •

All properties that included the behavior of the mode logic had to be made
more precise. At this point, the requirements were precise enough to have a
unique formal representation as (temporal) logic formulas over signals of the
TCM.

  * •

We defined five sub-requirements (e.g., A1 and A2 in Fig. 3) for properties
that were verified compositionally.

All requirements except the six for the mode logic and two requirements on
operational limits required making assumptions about the physics of the
aircraft.

##  4 Verification of the TCM

This section describes our verification efforts for the TCM: how we handle
Simulink models, and how we encode and verify properties. The complete TCM
benchmark (Simulink models and Lustre code) can be found in [1].

###  4.1 Handling Simulink Models

To apply SMT-based model checking, we compile the TCM Simulink model into the
synchronous dataflow language Lustre [6, 16]. In the following, we briefly
introduce Lustre and then describe the compilation process and how safety
properties are encoded and verified.

#### Lustre

Synchronous languages are a class of languages proposed for the design of
reactive systems (i.e., systems that maintain a permanent interaction with
their physical environment). Such languages are based on the theory of
synchronous time, in which the system and its environment are considered to
both view time with some “abstract” universal clock. Lustre combines each data
stream with an associated clock as a means to discretize time. The overall
system is considered to have a universal clock that represents the smallest
time span the system is able to distinguish, together with additional,
coarser-grained, user-defined clocks. Therefore the overall system may have
different subsections that react to inputs at different frequencies. At each
clock tick, the system is considered to evaluate all streams, so all values
are considered stable for any actual time spent in the instant between ticks.
Lustre programs and subprograms are expressed in terms of Nodes. Nodes
directly model subsystems in a modular fashion, with an externally visible set
of inputs and outputs. A node can be seen as a mapping of a finite set of
input streams (in the form of a tuple) to a finite set of output streams (also
expressed as a tuple).

#### Simulink to Lustre

In Matlab/Simulink from MathWorks©111http://www.mathworks.com/, dynamic
systems are modeled as block diagrams. Simulink uses dataflow-oriented block
diagram notation which consists of blocks and lines. Blocks represent either
some kind of functionality, like mathematical or logical functions, or they
are used for structuring the model in terms of subsystem blocks, port blocks,
bus blocks etc. Every block is defined by its type and its block parameters.

We have developed a tool called GAL222https://bitbucket.org/lememta/gal
(GeneAuto for Lustre) based on the _GeneAuto_ 333http://www.geneauto.org tool
set. The latter is a tool for the automatic code generation of Simulink models
to generate C, VHDL and Ada code. Although the development of the generator
was of primary interest, the GeneAuto project also put an emphasis on the
qualification of the toolchain [29] by providing traceability information all
along the code generation process.

Figure 2: Simulink model of the Altitude Controller subsystem with a safety
property encoded as synchronous observer.

GAL can only translate a subset of Simulink blocks. This subset can be
characterized as a collection of the most basic discrete-time blocks in
Simulink. In a typical controller model built by control engineers, one is
likely to encounter additional blocks such as the transfer function block, the
saturation blocks, the dead-zone block, and the integrator block. These blocks
have to first be transformed into equivalent Simulink models that GAL can
handle. We have developed an automated pre-processor for Simulink models, that
transforms an arbitrary Simulink model into a model “digestible” by GAL. The
pre-processor also generates an equivalence check between the original and
transformed models, modulo finite-precision arithmetic and discretization. The
check can be performed using standard simulation techniques supported by
Matlab. Note that our pre-processor provides automated support for
transformations that are standard among aerospace engineers when using the
MathWorks© Simulink Coder444http://mathworks.com/products/simulink-coder/.

#### Encoding safety properties

An extensively used technique to define expected behavior is synchronous
observers [17]. Synchronous observers provide an alternative to temporal
logics for specifying safety properties; the benefit of observers is that they
express properties in the same notation as the system model [26]. Observers
are typically used for simulation and testing purposes. A synchronous observer
is a wrapper used to test observable properties of a node $N$ with minimal
modification to the node itself; it returns an error signal if the property
does not hold. The task of checking the property is thus reduced to simply
checking if the stream is constantly true.

Synchronous observers are expressed in Simulink using a masked subsystem
block. A subsystem block is a container for a set of blocks. Masking a block
means extending it with some additional parameters. An example of such
synchronous observers expressed in Simulink is given in Figure 2. The red
blocks encode the safety property expressed as synchronous observer. GAL
translates these blocks as a property annotation to be proven in Lustre.

Specifically, a Lustre observer is a node taking as input all the flows
relevant to the safety property to be specified, and computing a Boolean flow
(e.g., “Obs” in Fig. 3) which is true as long as the observed flow satisfies
the property. We have used PKind [22] to prove the safety properties of Lustre
programs. PKind is a parallel k-induction-based model checker [22], which
includes automated invariant generation based on templates [21] and abstract
interpretation [13].

###  4.2 Safety verification results

Table 1 summarizes the verification results of the 20 safety properties on the
TCM model. At the beginning of the verification process, we discovered several
modeling errors within the TCM, which led to the falsification of some
properties:

  1. 1.

Some components produced output when disabled (e.g., the altitude controller).
This happened because the TCM model given to us was incomplete: the mode logic
was not implemented completely (which also affected the mode logic
properties). We remedied these problems by incorporating the necessary mode
logic into the model.

  2. 2.

Manual inputs from the pilot did not override the outputs of the autopilot for
all three axes. This was again due to an incompleteness in the TCM model. We
added a Simulink block (before the final output of the autopilot) to reflect
the fact a pilot has the ability to override the autopilot output.

  3. 3.

Some inputs were not variables but appeared as fixed constant values in the
model (e.g., the bank angle limit of G-240). This was simply a modeling error,
and it was easily corrected by modifying the appropriate variables.

  4. 4.

G-180 had to be refined to resolve a conflict with G-210 and the implicit
assumption that only the FPA control or the altitude control can be active at
any moment in time.

The results obtained for verification after the above changes are described in
Table 1. The properties are colored according to the verification technique
used. Gray properties (7-17) are the ones proved via a compositional argument.
Green properties (1-6) are the ones proved with a direct (non-compositional)
proof technique. Red properties (17-20) are the ones that could not be proven
on this specific model of the TCM. The latter properties applied to the B737
vehicle (see Section 3), however they referenced functionalities not
implemented in the TCM. It took an average of 2secs for PKind to verify the
green properties.

All properties were first attempted directly, without a compositional
argument. PKind was unable to verify the gray properties, despite a very high
timeout setting (5 hours). Since $k$-induction is sound but not complete, this
result has two possible interpretations: these properties are $k$-inductive
for an extremely high $k$ or, more plausibly, they are not $k$-inductive for
any $k$. We have also tried other Lustre verification tools based on different
verification techniques. Specifically, we used
Kind-2555http://kind2-mc.github.io/kind2/ and
Zustre666www.bitbucket.org/lememta/zustre. Kind-2 is a complete re-
implementation of PKind that also adds a verification engine based on IC3 [4],
while Zustre is a tool based on the generalized property-oriented reachability
implemented in Z3 [18]. Both tools were not able to prove the gray properties.

In a second step, we decomposed the properties either in terms of component-
level properties (e.g. G-120, G-130, G-140, G-150), or in “simpler” properties
to deal with (e.g. G-250, G-110). We now give details of the compositional
analysis of G-120 (see Section 3 for the natural language description). This
property involves 3 components of the TCM longitudinal control system: the
Mode Logic, the Altitude controller and the Flight Path Angle (FPA)
controller. In order to prove G-120, we decomposed the property into 4
component-level properties: G-180 and A1 for the Mode Logic; A2 for the
Altitude controller and FPA1 for the FPA controller. After proving the
component level properties, one still needs to make a formal compositional
argument that these properties imply the system-level property G-120. The
latter argument is captured in Figure 3. The upper box shows the Lustre nodes
of the various components involved in G-120. Each component comes with its own
guarantees. Such guarantees are used as assumptions in proving G-120
(described in the lower box).

⬇

node AutoPilot(HeadMode, AilStick, ElevStick, AltMode: real;

FPAMode, ATMode, AltCmd, Altitude, CAS,CASCmdMCP: real;)

returns (HeadEng, AltEng, FPAEng, ATEng: bool; CASCmd: real);

let

– G-180

assert (FPAMode = 0.0) or (not (AilStick = 0.0)) or

(not (ElevStick = 0.0)) or (FPAEng = true);

– A1

assert (not (AltMode= 0.0)) or (AltEng = false);

tel

node AltitudeControl (AltEng: bool; AltCmd, Alt: real;

GsKts, Hdot, HdotChgRate: real)

returns (AltGammaCmd: real ) ;

let

– A2

assert (AltEng = true) or (AltGammaCmd = 0.0);

tel

node FPAControl(Engage: bool; AltGammaCmd, Gamma: real;

ThetaDeg, VT: real)

returns (PitchCmd, PrePitchCmd: real);

let

– FPA1

assert true -> (Engage = false) or (AltGammaCmd = Gamma)

or((AltGammaCmd > Gamma) and (PitchCmd > pre(PrePitchCmd)))

or((AltGammaCmd < Gamma) and (PitchCmd < pre(PrePitchCmd)));

tel

⬇

node G-120 (HeadMode, AilStick, ElevStick, AltMode: real;

FPAMode, ATMode, AltCmd, Altitude, CAS: real;

CASCmdMCP, Gskts, Hdot, HDotChgRate, GammaCmd: real;

Gamma, ThetaDeg, VT: real)

returns (Obs: bool);

var

AltGammaCmd, FPain, TAlt, TFpa : real;

HeadEng, AltEng, FPAEng, ATEng : bool;

CasCmd, PitchCmd, PrePitch : real;

let

HeadEng,…,CasCmd = AutoPilot(HeadMode,…,CASCmdMCP);

AltGammaCmd = AltitudeControl(AltEng,…,HDotChgRate);

PitchCmd, PrePitch = FPAControl(FPAEng,…,VT);

assert FPain = (AltGammaCmd + GammaCmd);

assert (AltMode = 0.0);

assert (not (FPAMode = 0.0));

assert (ElevStick = 0.0);

assert (AilStick = 0.0);

assert (GammaCmd > 1.0 and GammaCmd < 10.0);

Obs = true -> (GammaCmd = Gamma)

or ((GammaCmd > Gamma) and (PitchCmd > pre(prePitch)))

or ((GammaCmd < Gamma) and (PitchCmd < pre(PrePitch)));

–!PROPERTY: obs = true;

tel

Figure 3: Compositional argument for G-120 property.

##  5 Lessons Learned

In this section, we summarize our experience and lessons learned from the
application of formal methods to the TCM case study. Some findings confirmed
our expectations: for example, we anticipated the fact that we would need to
consult with domain experts both for requirements elicitation, and for
assumption generation. We had, however, underestimated the extent to which
this would be required. Others surprised us: it took an extremely long time to
identify a case study that is both representative of flight-critical systems,
and publicly available; the requirements elicitation and formalization phases
were also much more involved than we expected. Regarding how far we would be
able to go with this case study, we had no expectations to begin with, since
we knew nothing about the system that we analyzed when we started our work.

1\. Case studies are hard to find. It is difficult to obtain real case studies
that are not proprietary and that can be shared outside an institution. The
TCM is available for General Purpose Release from NASA Langley, with case
number of LAR-18322-1. The process for obtaining the code is detailed in the
latest NASA Software
Catalog777http://technology.nasa.gov/NASA\\_Software\\_Catalog\\_2014.pdf. The
TCM benchmark used in this paper can be found in [1].

2\. System description must be massaged. Despite the progress in automating
verification techniques, a huge amount of effort still needs to be placed in
such a task. First of all, it is rare that verification tools are able to
directly handle all the features of the languages in which systems are
expressed. Despite the race towards keeping verification tools up-to-date,
modeling or programming languages are typically a step ahead. There is
therefore always an initial step involved, where the system description is
massaged to be handled by the targeted verification tools.

3\. Requirements elicitation. Requirements of flight-critical systems are
often hard to identify [15, 14]. Even when requirements are available, they
are very often written in natural language and need to be translated into a
clear notation with unambiguous semantics. We had seriously under-estimated
the effort that would be required in coming up with requirements that we could
verify. We involved collaborators from NASA’s Armstrong Flight Research
Center, and still required a lot of additional effort to bridge high-level
requirements with verifiable, component-level ones. It is the first time that
we had to tackle requirements starting from FARs, and we appreciate the
complexity of certification tasks, and the use of safety cases to organize
them [8] (we used safety case tools to organize our requirements, but cannot
present this work here due to limited space).

4\. Incomplete requirements and assumptions. Our case study confirmed that
requirements are often incomplete or even wrong. Developers often make
assumptions about the physics or the environment of a system that are not
explicitly expressed, and without which requirements do not hold. The
capability to analyze requirements with automated tools is invaluable in
identifying such problems with requirements. For example, the analysis we
performed for property G-120 (see Section 4) revealed the fact that an
implicit assumption needed to be formalized and become part of the
requirements.

5\. Scalability. The amount of progress in automating verification has been
substantial. However, there will always be properties and systems on which the
verification does not scale, unless more sophisticated compositional
approaches are introduced to break the problem into smaller, more manageable,
tasks. Compositional verification was needed to address 6 of the 20
requirements that we studied (gray properties in Figure 1). Decomposing
requirements was a non-trivial, manual task. Although our tools for this type
of system do not directly support automated assumption generation and
compositional verification yet, we believe that such techniques could
facilitate the application of our approaches [19].

6\. Domain expertise. Our case study proved that we are not yet at the point
where formal verification can occur in the absence of domain expertise. All of
the activities in observations 2, 3, and 4 above, required extensive
discussion with a domain expert that was part of the team assembled for this
study.

7\. Verification tools. Automated formal verification tools have made
tremendous progress in recent years. They are able to cope with the growing
complexities of systems. In our work, this was the key enabler in carrying out
the safety analysis of TCM. Specifically, SMT-based model checking was quite
effective in discharging the safety properties. In certain components of the
TCM we had to deal with nonlinear arithmetic operations (e.g., trigonometric
functions). While nonlinear arithmetic operations are not fully supported in
current SMT solvers, we were able to cope with that by using uninterpreted
functions 888The idea is to substitute nonlinear functions with uninterpreted
one. This could lead to non-feasible counterexamples. In this case we add
additional constraints to the uninterpreted function in order to eliminate
such counterexamples.. Handling nonlinear arithmetic operations is essential
for the verification of flight-critical systems; it is therefore desirable to
develop tools that can robustly handle these features.

9\. Level of effort. Verifying the TCM required approximately three person
months (see Table 2), with the involvement of verification and aeronautics
experts. Most of our time was spent eliciting the properties, formalizing
them, and creating the assumptions about the physical environment that we need
for verification. We also spent several weeks working towards utilizing our
verification within a future safety case effort. The actual verification
process was automated, and required the least time.

Table 2: TCM verification: Approximate level-of-effort Effort | Person Months
---|---
Implementation (Tools) | 0.50
Preparation of Models | 0.50
Property Elicitation | 0.25
Formalization (Relation to model) | 0.50
Physical Assumptions | 0.50
Compositional Arguments | 0.25
Verification | 0.20
Safety Case Generation | 0.50

##  6 Related Work

Our aim in designing this study was to make it as realistic, independent, and
shareable as possible, and hence we targeted a system that is representative
of flight-critical systems, that was developed outside our group, and that is
available to the research community. It is hard to find realistic studies in
the research community that are not proprietary and that can therefore be used
as benchmarks.

In recent work, we have applied probabilistic verification and synthesis
techniques to analyze the ACAS X onboard collision avoidance system [30].
Moreover, we have developed a testing infrastructure for the automated
analysis of the AutoResolver air-traffic control system, aimed at the
prediction and resolution of aircraft loss of separation [15]. A previous
large study performed by our group aimed at comparing model checking, static
analysis, runtime analysis and testing, through their application for finding
bugs in the Executive component of an autonomous robot developed at NASA Ames
[5]. These systems are not publicly available.

Several studies related to the verification of flight-critical systems have
been performed by Rockwell Collins. In [24], the authors report on the use of
their automated framework to verify Simulink and Stateflow designs of three
aeronautics components: the ADGS-2100 Window Manager, for ensuring that data
from different applications is routed to the correct aircraft display panel;
two components of the operational flight program of an unmanned aerial vehicle
developed by Lockheed Martin Aerospace, one involving redundancy management,
and the other one in charge of generating actuator commands for the aircraft’s
six control surfaces. These studies confirmed the applicability and benefits
of formal verification techniques in the design of flight-critical systems.
The commercial components that were targeted are not publicly available.

More recently, Rockwell Collins has developed compositional techniques for
scalable verification of architectural models expressed in the AADL language
[7]. Funded through a NASA NRA, they have recently applied compositional
verification to TCM architectural models, in the context of the same program
milestone that drove this work999The case studies performed in this work can
be found in www.github.co/smaccm.

The Astrée static analyzer has been used to prove the absence of runtime
errors from two Airbus components implemented in C [28, 2], as well as from a
C version of the automatic docking software of the Jules Vernes Automated
Transfer Vehicle (ATV) enabling ESA to transport payloads to the International
Space Station [3]. Galdino et al [11] use the PVS theorem prover to formally
verify an air-traffic control resolution and recovery algorithm. In the domain
of hybrid system verification, Platzer and Clarke [25] have applied the
KeYmaera verification tool to prove properties of curved flight collision
avoidance maneuvers. Esteve et al. have applied a probabilistic model checker
to determine properties of an early design spacecraft model for the European
Space Agency [9].

For an extensive study of success stories related to the application of formal
verification in practice, we refer the reader to the following technical
report by Garavel and Graf [12]. Note that, in this paper, we focus on case
studies related to avionics; several other studies of safety-critical systems
have been performed, for example in the contexts of medical devices and of the
automotive industry.

##  7 Conclusion

To summarize, we demonstrated a verification approach for the TCM controls
system: a publicly available, realistic and complex flight-critical system of
moderate size. This study required a significant amount of effort from a team
made up of both verification and domain experts. Compositional verification
was required to prove some of the safety properties of the system. The only
safety properties we did not prove in this study were those in which the
desired functionality had not actually been modeled. Our experience highlights
the promise of compositional verification in the certification of flight-
critical systems. In practice, we saw that the most significant part of our
effort was in defining and formalizing the appropriate properties, starting
from high-level FARs requirement, all the way down to properties of the target
system. We hope that our experience will be useful to other researchers and
will encourage them to apply and evaluate alternative techniques on the TCM
case study. To this aim, we have made all the artefacts of our benchmark
publicly available.

## References

  * [1] TCM models and lustre code.  Avaiable on https://bitbucket.org/lememta/tcm_comp_verification/wiki/Home.
  * [2] B. Blanchet, P. Cousot, R. Cousot, J. Feret, L. Mauborgne, A. Miné, D. Monniaux, and X. Rival.  Design and implementation of a special-purpose static program analyzer for safety-critical real-time embedded software, invited chapter.  In T. Mogensen, D. Schmidt, and I. Sudborough, editors, The Essence of Computation: Complexity, Analysis, Transformation. Essays Dedicated to Neil D. Jones, LNCS 2566, pages 85–108. Springer-Verlag, Oct. 2002.
  * [3] O. Bouissou, E. Conquet, P. Cousot, R. Cousot, J. Feret, K. Ghorbal, E. Goubault, D. Lesens, L. Mauborgne, A. Miné, S. Putot, X. Rival, and M. Turin.  Space software validation using abstract interpretation.  In Proc. of the Int. Space System Engineering Conf., Data Systems in Aerospace (DASIA 2009), volume SP-669, pages 1–7, Istambul, Turkey, May 2009. ESA.
  * [4] A. R. Bradley.  SAT-based model checking without unrolling.  In Proceedings of the 12th International Conference on Verification, Model Checking, and Abstract Interpretation, VMCAI’11, pages 70–87. Springer-Verlag, 2011.
  * [5] G. P. Brat, D. Drusinsky, D. Giannakopoulou, A. Goldberg, K. Havelund, M. R. Lowry, C. S. Pasareanu, A. Venet, W. Visser, and R. Washington.  Experimental evaluation of verification and validation tools on Martian rover software.  Formal Methods in System Design, 25(2-3):167–198, 2004.
  * [6] P. Caspi, D. Pilaud, N. Halbwachs, and J. A. Plaice.  Lustre: a declarative language for real-time programming.  In Proceedings of the 14th ACM SIGACT-SIGPLAN symposium on Principles of programming languages, POPL ’87, pages 178–188. ACM, 1987.
  * [7] D. D. Cofer, A. Gacek, S. P. Miller, M. W. Whalen, B. LaValley, and L. Sha.  Compositional verification of architectural models.  In NASA Formal Methods - 4th International Symposium, NFM 2012, Norfolk,VA, USA, April 3-5, 2012. Proceedings, pages 126–140, 2012.
  * [8] E. Denney, G. Pai, and J. Pohl.  AdvoCATE: An assurance case automation toolset.  In SAFECOMP Workshops, pages 8–21, 2012.
  * [9] M. Esteve, J. Katoen, V. Y. Nguyen, B. Postma, and Y. Yushtein.  Formal correctness, safety, dependability, and performance analysis of a satellite.  In 34th International Conference on Software Engineering, ICSE 2012, June 2-9, 2012, Zurich, Switzerland, pages 1022–1031, 2012.
  * [10] Federal Aviation Administration.  Electronic code of federal regulations.
  * [11] A. L. Galdino, C. M. noz, and M. Ayala-Rincón.  Formal verification of an optimal air traffic conflict resolution and recovery algorithm.  In Logic, Language, Information and Computation, 14th International Workshop, WoLLIC 2007, Rio de Janeiro, Brazil, July 2-5, 2007, Proceedings, pages 177–188, 2007.
  * [12] H. Garavel and S. Graf.  Formal methods for safe and secure computer systems.  Technical Report BSI-Study 875, Bundesamt fuer Sicherheit in Informationstechnik, December 2013.
  * [13] P.-L. Garoche, T. Kahsai, and C. Tinelli.  Incremental invariant generation using logic-based automatic abstract transformers.  In NFM, pages 139–154, 2013.
  * [14] D. Giannakopoulou, D. H. Bushnell, J. Schumann, H. Erzberger, and K. Heere.  Formal testing for separation assurance.  Ann. Math. Artif. Intell., 63(1):5–30, 2011.
  * [15] D. Giannakopoulou, F. Howar, M. Isberner, T. Lauderdale, Z. Rakamaric, and V. Raman.  Taming test inputs for separation assurance.  In 19th IEEE/ACM International Conference on Automated Software Engineering (ASE 2014), 2014.
  * [16] N. Halbwachs, P. Caspi, P. Raymond, and D. Pilaud.  The synchronous dataflow programming language Lustre.  In Proceedings of the IEEE, pages 1305–1320, 1991.
  * [17] N. Halbwachs, F. Lagnier, and P. Raymond.  Synchronous observers and the verification of reactive systems.  In AMAST, pages 83–96, 1993.
  * [18] K. Hoder and N. Bjørner.  Generalized property directed reachability.  In A. Cimatti and R. Sebastiani, editors, Theory and Applications of Satisfiability Testing – SAT 2012, volume 7317 of LNCS, pages 157–171. 2012.
  * [19] F. Howar, D. Giannakopoulou, and Z. Rakamaric.  Hybrid learning: interface generation through static, dynamic, and symbolic analysis.  In International Symposium on Software Testing and Analysis, ISSTA ’13, Lugano, Switzerland, July 15-20, 2013, pages 268–279, 2013.
  * [20] R. M. Hueschen.  Development of the Transport Class Model (TCM) aircraft simulation from a sub-scale Generic Transport Model (GTM) simulation.  Technical report, NASA, Langley Research Center, Hampton, VA, August 2011\.
  * [21] T. Kahsai, Y. Ge, and C. Tinelli.  Instantiation-based invariant discovery.  In NFM, volume 6617 of LNCS, pages 192–207, 2011.
  * [22] T. Kahsai and C. Tinelli.  PKIND: a parallel $k$-induction based model checker.  In PDMC, volume 72 of EPTCS, pages 55–62, 2011.
  * [23] A. Lyapunov.  General problem of the stability of motion.  PhD thesis, Univ. Kharkov, 1892.
  * [24] S. P. Miller, M. W. Whalen, and D. D. Cofer.  Software model checking takes off.  Commun. ACM, 53(2):58–64, 2010.
  * [25] A. Platzer and E. M. Clarke.  Formal verification of curved flight collision avoidance maneuvers: A case study.  In A. Cavalcanti and D. Dams, editors, FM 2009: Formal Methods, 16th International Symposium on Formal Methods, volume 5850 of LNCS, pages 547–562. Springer, 2009.
  * [26] J. Rushby.  The versatile synchronous observer.  In 15th Brazilian conference on Formal Methods: foundations and applications, SBMF’12, pages 1–1, Berlin, Heidelberg, 2012. Springer-Verlag.
  * [27] SmartCockpit.  B737 automatic flight systems summary.
  * [28] J. Souyris and D. Delmas.  Experimental assessment of Astrée on safety-critical avionics software.  In SAFECOMP, pages 479–490, 2007.
  * [29] A. Toom, N. Izerrouken, T. Naks, M. Pantel, and O. Ssi-Yan-Kai.  Towards reliable code generation with an open tool: Evolutions of the Gene-Auto toolset.  In ERTS, http://www.sia.fr, 2010. Société des Ingénieurs de l’Automobile.
  * [30] C. von Essen and D. Giannakopoulou.  Analyzing the next generation airborne collision avoidance system.  In Tools and Algorithms for the Construction and Analysis of Systems - 20th International Conference, TACAS 2014, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2014, Grenoble, France, April 5-13, 2014. Proceedings, pages 620–635, 2014.


Breast Tumor Cellularity Assessment using Deep Neural Networks

Alexander Rakhlin

Neuromation OU

Tallinn, 10111 Estonia

Aleksei Tiulpin

University of Oulu

Oulu 90220, Finland

Alexey A. Shvets

Boston, MA 02142, USA

Alexandr A. Kalinin

University of Michigan

Ann Arbor, MI 48109, USA

Vladimir I. Iglovikov

San Francisco, CA 94107, USA

Sergey Nikolenko

Neuromation OU, Estonia,

Steklov Mathematical Institute
at St. Petersburg, Russia

Breast cancer is one of the main causes of death worldwide. Histopathological cellularity assessment of residual tumors in post-surgical tissues is used to analyze a tumor's response to a therapy. Correct cellularity assessment increases the chances of getting an appropriate treatment and facilitates the patient's survival. In current clinical practice, tumor cellularity is manually estimated by pathologists; this process is tedious and prone to errors or low agreement rates between assessors. In this work, we evaluated three strong novel Deep Learning-based approaches for automatic assessment of tumor cellularity from post-treated breast surgical specimens stained with hematoxylin and eosin. We validated the proposed methods on the BreastPathQ SPIE challenge dataset that consisted of 2395 image patches selected from whole slide images acquired from 64 patients. Compared to expert pathologist scoring, our best performing method yielded the Cohens kappa coefficient of $0.69$ (vs. $0.42$ previously known in literature) and the intra-class correlation coefficient of $0.89$ (vs. $0.83$). Our results suggest that Deep Learning-based methods have a significant potential to alleviate the burden on pathologists, enhance the diagnostic workflow, and, thereby, facilitate better clinical outcomes in breast cancer treatment.

§ INTRODUCTION

figureGeneric description of the methods developed and evaluated in this study. Our first approach leverages segmentation model, feature extraction and gradient boosted trees. The second approach directly predicts the cellularity from the raw data. Finally, in our third setting, we combine the first and the second approach and used a deep convolutional neural network to predict the cellularity score from segmentation mask.

Breast cancer is one of the most common cancer types diagnosed in women in the United States and worldwide [36]. Biopsies and histological assessment allow pathologists to analyze microscopic structures of breast tissues and, in particular, assess the cancer's aggressiveness.

Multiple options are available to manage and monitor the breast cancer treatment based on the information provided from the tumor's response to it. In addition to the treatment effect on the tumor size, the therapy may also alter the tumor's cellularity [8]. During anticancer therapy, the size of the tumor may remain the same, but the overall cellularity may be drastically reduced [30]. As a result, it makes the residual tumor cellularity an important factor in assessing the response treatment.

Currently, tumor cellularity is manually assessed by pathologists from hematoxylin and eosin (H$\&$E)-stained slides  [30]. The costs of such estimation are high, the process is tedious and subjective, and the quality and reliability might be also be affected by high inter-observer variability even among senior pathologists. This potentially may affect prognostic power assessment in clinical trials [39]. The subjectivity in visual tissue assessment motivates the use of computer-aided methods to improve the diagnosis accuracy, reduce human error and increase inter-observer agreement and reproducibility [27, 11]. Automated analysis of the H&E slide using computer vision could provide immediate benefits to patient care. Recent success in Deep Learning (DL) [22, 34], and in particular the advances in convolutional neural networks (CNN), have recently shown high potential in this realm [9].

In this work, we evaluate three DL-methods to score the cellularity of the breast tissue from histopathological images. In particular, our first approach employs a weakly-supervised segmentation model with Resnet-34 [13] encoder and Feature Pyramid Network (FPN)[24] and a second-stage regression network that predicts the cellularity score using the predicted segmentation maps. Our second approach is also based on segmentation, however, instead of using the segmentation maps directly, we extract various features from them and use the gradient boosting trees (GBT) [28] to predict the cellularity score. Finally, we also evaluate using H&E image patches directly to predict the cellularity score.

§ RELATED WORK

CNNs have recently been successfully applied to many tasks in biomedical image analysis, often outperforming conventional machine learning methods [9, 41, 35]. As such, they have successfully been utilized for digital pathology image analysis and have demonstrated great potential for improving breast cancer diagnostics [38, 4, 5, 32].

Although there are not many studies focusing directly on automated quantitative cellularity assessment, it has been shown that this task can be solved by first segmenting malignant cells and then computing the tumor's area [29]. Many efforts have been devoted to developing supervised and unsupervised methods for automated cell and nuclear segmentation and detection [44, 21]. Supervised segmentation models have superior performance but require hand-labeled nuclear mask annotations [44]. In these approaches, segmented nuclear bodies are used to extract features that are typically inspired by visual markers recognized by pathologists. Commonly used features describe morphology, texture, and spatial relationships among cell nuclei in tissue [29, 19].

The conventional approach most relevant to our work is by Peikari et al. [29] who proposed an automated cellularity assessment protocol. First, they used smaller patches, or regions of interest (RoI), extracted from whole slide images to segment all present cell nuclei. Then they extracted a number of predefined features from segmented nuclei and used support vector machines to distinguish lymphocytes and normal epithelial nuclei from malignant ones. Cellularity estimation was done using distinguished malignant epithelial figures for every RoI.

Alternatively, segmentation-free methods that directly estimate cellularity from histopathology imaging data and nuclei locations annotated by human observers were also shown promising. In particular, Veta et al. [43] proposed a deep learning-based method that leverages an information from a tumor's cells nuclei locations (centroids) and predicts the areas of individual nuclei and mean nuclear area without the intermediate step of nuclei segmentation. In particular, this approach was based on a 10-layer deep neural network predicting nuclear areas quantized into 20 histogram bins. The results showed that predicted measurements had substantial agreement with manual measurements, which suggests that it is possible to compute the areas directly from imaging data, without the intermediate step of nuclei segmentation. This is in spirit similar to one of our approaches, but we do not directly compare our methods to Veta et al. since we use different datasets and performance metrics.

figures/fpn-net-crop.pdf18pdfcrop figures/fpn-net.pdf

Encoder-decoder segmentation network architecture with Resnet-34 encoder and feature pyramid network decoder. Spatial Dropout 2D is added after multi-layer concatenation.

Recent works by Akbar et al. [3, 2] have compared the conventional approach based on segmentation and feature extraction and direct applications of deep CNNs to image patches in both regression and classification settings. Overall, they showed that the DL-based approach outperformed hand-crafted features in both accuracy and intra-class correlation (ICC) with expert pathologist annotations. Specifically, their best result was achieved by using a pre-trained Inception [40] model that reached ICC of $0.83$ and $0.81$ with two expert pathologists. In this study we evaluate even wider range of DL-based approaches, including segmentation-based and segmentation-free, in both regression and classification settings. We provide appropriate performance comparisons with previously reported results[It is worth noting that the official BrestPathQ challenge results have been reported only as a score distribution. Each team know  their  own  results  only, and ours belong to the  right end of the distribution, but, unfortunately,  we are not able to provide a comparison of our approach  with other participants. <http://spiechallenges.cloudapp.net/competitions/14#learn_the_details>]. All the methods developed in this study are fully automatic and do not require any involvement of the human annotators at the test time.

§ METHODS

In this study, we propose and evaluate three different methods. The first two methods are based on the nuclei segmentation and the third method leverages the raw image without preceding segmentation step. Graphical illustration of our approach is presented in Figure <ref>

§.§ Segmentation

Network Architecture. Most modern segmentation architectures inherit the encoder-decoder architecture similar to U-Net [33], where convolutional layers in the contracting branch (encoder) are followed by an upsampling branch that brings segmentation back to the original image size (decoder). In addition, skip connections are used between contracting and upsampling modules to help the localization information propagate through the complex multilayer structure and eventually improve segmentation accuracy [33]. U-Net and architectures inspired by this idea have produced state of the art results in various segmentation problems, and many improvements for the architecture and its training protocols have recently been proposed. In particular, Iglovikov et al. [15] used batch normalization [17] and exponential linear unit (ELU) as the primary activation function and an ImageNet pre-trained VGG-11 network [37] as an encoder. Liu et al. [25] proposed an hourglass-shaped network (HSN) with residual connections, which is also very similar to the U-Net architecture. Rakhlin et al. [31] used the Resnet-34 network [14] as the encoder and the Lovász-Softmax loss function [6] along with Stochastic Weight Averaging (SWA) [18] for training.

In our proposed architecture, the segmentation module also inherits the U-Net architecture. The contracting branch (encoder) of our model is based on the Resnet-34 [14] network architecture where we have introduced several useful modifications. In particular, we have replaced ReLU activations with ELU that does not saturate gradients and keeps the output close to zero mean and have changed order of batch normalization [17] and activation layers. In Section <ref> we compare encoders initialized with random He's initialization [12] and pretrained on ImageNet.

Light micrograph of a histologic specimen of breast tissue stained with hematoxylin and eosin (top). The bottom row shows nuclei segmentation masks synthesized from weak labels: Malignant — red, Normal — green, Lymphocyte — blue.

To address the limited size of the BreastPathQ Cancer Cellularity Challenge dataset, we utilized two regularization techniques: (1) data augmentation and (2) spatial 2D dropout incorporated into the upsampling branch [42].

The upsampling branch is implemented as a Feature Pyramid Network (FPN) [23], reconstructing high-level semantic feature maps at four scales simultaneously. We implement a feature pyramid block as a convolutional layer with 64 activation maps followed by upsampling to the original resolution with upsampling rate of 8, 4, 2, or 1 depending on the feature map depth (see Fig. <ref>). In Section <ref>, we compare the performance of standard and FPN decoders.

We concatenate upsampled maps into a single layer of $64\times4=256$ maps and add after it a spatial 2D dropout layer, which acts as a regularizer and prevents coadaptation of the network weights, but unlike conventional dropout it drops out not individual neurons but rather entire activation maps. Throughout the work, we use dropout rate $0.5$, randomly dropping $128$ out of $256$ activation maps.

Finally, the output of the model is a 4-channel sigmoid layer that assigns every pixel with four values from $0$ to $1$ that represent the probabilities of belonging to the Normal, Lymphocyte, Malignant, and Background classes.

Loss functions. Binary cross entropy (BCE), while convenient for training, does not directly translate into Jaccard index, the metric commonly used to evaluate segmentation accuracy. Hence, as the loss function we use

\begin{equation}
L^c(w)=(1-\alpha) \mathrm{BCE}^c(w) - \alpha J^c(w),
\label{loss}
\end{equation}

a weighted sum of $\mathrm{BCE}$ and the soft Jaccard loss for class $c$ [15, 31, 16]. In this work, we set $\alpha=0.15$, a value found via cross-validation. The soft Jaccard loss is defined as

\begin{equation}
\label{jaccard_loss}
\end{equation}

where $w$ are network parameters, $y^c_i$ is the binary label for pixel $i$ and class $c$, $\hat{y}^c_i$ is the predicted probability of $c$ for pixel $i$, and $N$ is the total number of pixels.
The total loss function is a weighted sum of class losses:

\begin{equation}
L(w)=\frac{1}{V}\sum\limits_{c=1}^4L^c(w)v^c, V=\sum\limits_{c=1}^4v^c,
\label{total_loss}
\end{equation}

where $v^c$ is a loss weight for class $c$. In this work we weigh Normal, Lymphocyte, and Background as $1$ and Malignant, the class of primary importance in our problem, as $4$.

§.§ Cellularity estimation from segmented cells
In this subsection, we describe the method for cellularity assessment that leverages the output of the trained segmentation network Fig. <ref>. We feed the segmented output into a Resnet-34 CNN model. The model automatically learns deep features from the 4-channel segmentation input and regresses it onto continuous cellularity score using continuous regression loss ($L_2)$. In this approach, the segmentation model acts as a filter the aim of which is to extract only the information about the cell morphology. We hypothesized that this structured approach makes our method similar to methods employed by expert pathologists, makes it transparent and less sensitive to data acquisition settings.

§.§ Feature extraction-based cellularity estimation

The second type of model is Gradient Boosted Trees (GBT) [20] in regression mode ($L_2$ loss). The general idea and handcrafted features are borrowed from the second place solution for 2017 Kaggle contest for Sea Lion Population Count in aerial imaginary [26]. The authors would like to thank Konstantin Lopuhin for valuable discussion they had while incorporating his method. In this study, GBT operates on a vector of hand-crafted features extracted from nuclei segmentation maps, including:

* activations and their areas aggregated over segmentation maps with different thresholds; for every segmentation map in Normal, Lymphocyte, Malignant and for 7 thresholds $0.02$, $0.04$, $0.08$, $0.16$, $0.24$, $0.32$, $0.5$, we obtain 2 values: total area above threshold and total activation above threshold (see Fig. <ref> for an illustration);

* using the Laplacian of Gaussian (LoG) method as implemented in the OpenCV library [7], we find blobs in segmentation maps at 6 thresholds: $0.02$, $0.04$, $0.08$, $0.16$, $0.24$, $0.5$; for each threshold we find the number of blobs and total activation in blob centers (Fig. <ref>);

* total activation for every channel, computed as a sum of the activations at every pixel after sigmoid.

In total, we obtain $3 \times (7 \times 2 + 6 \times 2 + 1) = 81$ features to train the GBT model.

Segmentation results at thresholds (a) - (d) of Malignant channel superimposed with the original image. Masks generated after thresholding were used for feature extraction.

Nuclei blobs detected from the Malignant segmentation maps using the Laplacian of Gaussian method at thresholds (a)-(d). The blobs were used for feature extraction.

§.§ Cellularity estimation from the raw images

The third type of model is a deep convolutional network implemented in regression or classification settings ($L_2$ or categorical cross-entropy loss functions respectively). These models do not use intermediate segmentation and predict the cellularity score immediately from the microscopic image. For classification, we categorize cellularity into $101$ class using regular bins with thresholds $0.00,0.01,\ldots,1.00$. The idea of direct regression of an image into continuous value using CNN is not new. In particular, it was implemented in  [16] where the authors use CNN to predict bone age from radiograph.

§.§ Evaluation metrics

We assessed the results using several metrics. The main evaluation metric is the mean squared error (MSE) between the cellularity score obtained in our experiments and ground truth provided by an expert pathologist.

In order to make our results comparable with previous work, we also report Cohen's kappa coefficient agreement and the intra-class correlation coefficient (ICC) between expert and automated methods, similar to [29]. In all experiments, we find our results superior to our predecessors; however, the cellularity score itself in [29] is evaluated based on binning it into four categories of $0\mbox{--}25$%, $26\mbox{--}50$%, $51\mbox{--}75$%, and $76\mbox{--}100$%. Such $4$-class categorization is relatively coarse and, in our opinion, does not represent a suitable evaluation metric for continuous cellularity estimation that is our goal in this work.

§ EXPERIMENTS AND RESULTS

§.§ Data

The data used in this study had been acquired from the Sunnybrook Health Sciences Centre with funding from the Canadian Cancer Society and was made available for the BreastPathQ challenge sponsored by the SPIE, NCI/NIH, AAPM, and the Sunnybrook Research Institute [29].

every axis/.append style=
y label style=at=(axis description cs:0.1,.5),anchor=south,
legend style=
anchor=south east,
legend columns=1,
legend cell align=left,

every axis plot/.append style=
line width=1.pt,
mark size=1.25pt,
mark repeat=24,

In our experiments we used $2,395$ patches of $512\times512$ pixels in size, extracted from $96$ haematoxylin and eosin (H&E) stained whole slide images (WSI) acquired from 64 patients. Each patch in the training set has been assigned a tumor cellularity score by an expert pathologist. In Figure <ref>, we present a distribution of the cellularity scores in the dataset.

Besides the image data, we used the annotations ($X$ and $Y$ coordinates) to identify lymphocytes, malignant epithelial, and normal epithelial cell nuclei in the additional $153$ patches. Using these weak annotations, we generated the segmentation masks that were used in our experiments. Here, at each $XY$ location, we simply fit a blob of $15$ pixels in diameter. In Figure <ref> we present the generated masks for various classes.

xtick=0,0.2,0.4,0.6,0.8,1, xticklabels=0,0.2,0.4,0.6,0.8,1,
[ybar interval,fill=blue!80!black,opacity=.4] coordinates
(0,806) (0.07142857,193) (0.14285714,261) (0.21428571,100) (0.28571429,117) (0.35714286,171) (0.42857143,19) (0.5,164) (0.57142857,102) (0.64285714,109) (0.71428571,30) (0.78571429,83) (0.85714286,61) (0.92857143,178) (1,0);

Cellularity score distribution.

§.§ Weakly-supervised cell segmentation

In this study, characteristic features of the data present a serious challenge for developing segmentation models: 1) the data features no segmented cells, only their coordinates (that is why we use semi-supervised segmentation); 2) annotated nuclei are present only in 154 microscopic images, each containing 0-50 malignant cells. However, cell segmentation is not a distinct goal of this study. As mentioned in Section <ref> and in [8, 30, 10], cellularity within the tumor area is assessed by estimating the percentage area of the overall tumor bed comprised of invasive tumor cells.  Aggregated area of individual invasive cell areas serves as a proxy and does not represent the ultimate cellularity value. Cellularity is affected by cell density, localization, and tissue structure.  We use segmented cells essentially as an interpretable visualization of an invasive tumor within the tumor bed.

In our ablation studies, we evaluated our model in four different settings to find how different design choices influence the segmentation accuracy and generalization. Namely, we compared the model as described in Section <ref> with a standard U-Net decoder against an FPN decoder, and with the encoder initialized randomly against the encoder initialized with weights pretrained on ImageNet. In all settings, the model was trained for 150 epochs with the Adam optimizer and gradually decreasing learning rate from $10^{-4}$ to $10^{-5}$.

To obtain training patches, we downscaled the microscopy images $\times2$ times, randomly cropped a $256\times256$ area, and rescaled pixel values from $[0, 255]$ to $[-1, 1]$. As mentioned previously, segmentation targets were generated as 4-channel masks with round blobs, 15 pixels in diameter (the characteristic nucleus size), drawn in the nuclei centers. During training, we dynamically augmented images with vertical and horizontal flips, rotation, gamma, hue, and saturation utilizing the Albumentations library [1].

Segmentation results: the Jaccard index for different decoders and initializations.

Initialization     Standard decoder     FPN decoder
Random     0.35     0.47

ImageNet     0.50	    0.53

Cellularity MSE with 95% confidence intervals for the segmentation-based (first row) and for the end-to-end methods. Our results demonstrate the importance of ImageNet pre-training. C in the parentheses indicates classification, R – regression and S – segmentation.

2*Model     2cInitialization
Random     ImageNet
GBT     0.023 [0.019-0.026]     0.022 [0.019-0.026]

Resnet34 (SR)     0.013 [0.011-0.015]      0.013 [0.011-0.015]
ResNet34 (R)     0.015 [0.013-0.018]     0.011 [0.010-0.012]

ResNet50 (R)     0.025 [0.022-0.028]     0.011 [0.009-0.012]

Xception (R)     0.017 [0.015-0.020]     0.010 [0.009-0.012]

Xception (C)     0.012 [0.010-0.014]     0.010 [0.009-0.012]

Cellularity Kappa (4 class binning) and Intra-Class Correlation Coefficient (ICC) with 95% confidence intervals for the segmentation-based ($1^{st}$ and $2^{nd}$ rows) and for the methods predicting cellularity directly, without segmentation. All the models here utilize ImageNet pre-training. C in the parentheses indicates classification, R – regression and S – segmentation.

2*Model     2cMetric
Kappa     ICC
GBT     0.571 [0.520-0.622]     0.787 [0.744-0.823]
Resnet34 (SR)     0.658 [0.604-0.704]     0.865 [0.835-0.891]
ResNet34 (R)     0.649 [0.599-0.700]     0.868 [0.840-0.892]

ResNet50 (R)     0.652 [0.603-0.701]     0.867 [0.844-0.894]

Xception (R)     0.669 [0.616-0.713]     0.881 [0.853-0.904]

Xception (C)     0.689 [0.642-0.734]     0.883 [0.858-0.905]

Peikari et al. [29]     0.38-0.42     0.75 [0.71-0.79]

Akbar et al. [3]     —     0.83 [0.79-0.86]

Examples of the generated segmentation masks in the Malignant channel. Left to right: (a) original patch; (b) ground truth segmentation superimposed on the original image; (c) activation map; (d) nuclei blobs reconstructed from the activation map with the LoG method.

In the first series of experiments, we evaluated segmentation quality as an important intermediate metric for the evaluation of our methods. The segmentation performance as a function of the decoder and initialization is shown in Table <ref>. As we can see, the model with the feature pyramid decoder and encoder pretrained on ImageNet achieved significantly higher and more stable Jaccard index on the validation set than the alternatives. Figure <ref> shows an example of generated segmentation masks in the Malignant channel and nuclei blobs reconstructed with the Laplacian of Gaussian method.

[scaled y ticks=false,xmin=0,xmax=104,ymax=0.06,ymin=0,width=,height=.8,xlabel=Training epochs,ylabel=Cellularity score MSE,legend style=at=(1.0,1.0),anchor=north east,ytick=0,0.02,0.04,0.06,0.08,0.1, yticklabels=0,0.02,0.04,0.06,0.08,0.1,
y label style=at=(axis description cs:-0.1,.5)]
[smooth,mark phase=1,mark=*,color=blue!50!black] table[col sep=comma,x index=0,y index=3] mse.csv;
[smooth,mark phase=7,mark=o,color=red!50!black] table[col sep=comma,x index=0,y index=4] mse.csv;
[smooth,mark phase=13,mark=oplus,color=green!50!black] table[col sep=comma,x index=0,y index=1] mse.csv;
[smooth,mark phase=19,mark=triangle*,color=magenta] table[col sep=comma,x index=0,y index=2] mse.csv;
Training set, random init,Test set, random init,Training set, ImageNet init,Test set, ImageNet init

Cellularity MSE evolution during training.

§.§ Segmentation-based cellularity assessment

Prediction from the segmented cells. As mentioned previously, we used the output of the segmentation model as input for the cellularity regressor and then trained this cascade end-to-end. We froze the segmentation model and stack its 4-channel output with a randomly initialized Resnet-34 in the regression setting. We trained the regression part with cellularity targets and MSE loss until convergence. Then we unfroze segmentation weights and fine-tuned both modules in an end-to-end fashion, as a single model. We repeated this experiment with Resnet-34 pretrained on ImageNet. In the latter case, we excluded the background channel from segmentation output to comply with the vanilla Resnet-34 architecture that has a 3-channel input.

In these experiments, we found that after fine-tuning the accuracy of segmentation itself slightly decreases, while the accuracy of the overall cellularity scoring increases. This is in line with [10], which found that perfect segmentation of nuclei figures does not ensure better classification of malignant objects from breast cancer tissues. This finding suggests that the two branches of future work, tumor bed segmentation and cellularity assessment, are relatively independent.

Feature extraction-based method. In this series of experiments, we extracted the 81 features from segmentation masks as discussed in Section <ref> and trained the LightGBM [20] regression model with mean squared error (MSE) objective. The model was trained for 600 epochs with learning rate $0.01$. The maximum tree depth was set to 5; the number of leaves, to 8. These parameters have been selected through cross-validation.

We report LightGBM accuracy in Table <ref> and show the resulting feature importance on Fig. <ref>. Feature importance was calculated based on the total gain of the loss function from the splits formed according to this feature. As expected, all highly important features come from the Malignant channel. The most important feature is the total activation above $0.5$ threshold, and the second and third most important features are the activations above $0.32$ and $0.24$ thresholds, as expected since activations at different thesholds are highly correlated, and the segmentation quality at threshold $0.5$ was the best, so the feature based on this mask is a natural candidate for the most important feature. Activations at lower thresholds provide additional value, but a big part of the information that they contain has already been conveyed via the $0.5$ threshold feature. Interestingly, malignant cell count (detected at threshold $0.24$) is only the $9$th feature in order of importance.

§.§ Direct cellularity assessment from the raw images

In our final experiments, we evaluated several deep neural architectures that take the original microscopy images as input and output the cellularity score without intermediate segmentation. Similarly to previous experiments, we trained the models with random He initialization [12] or initialized them with weights pretrained on ImageNet. In all cases, ImageNet initialization was superior to random, and the overall accuracy was slightly better than for the models with intermediate segmentation. The Xception model implemented in a classification setup with random initialization performed slightly better than its counterparts (MSE 0.012 vs. 0.017-0.025). Although small, this difference could possibly be attributed to the known regression-to-mean problem of the continuous regression with $L_2$ loss (e.g., see a well-explained example for a colourization application in [45]). The mean squared error of cellularity prediction as a function of the training epoch for different initializations is shown in Figure <ref>. All the performance evaluation metrics are presented in Table <ref> and Table <ref>.

every axis plot/.append style=
line width=.6pt,
mark size=1.75pt,
mark repeat=10,

xbar, xmode=log, xtick=10,100,1000,10000,100000, xticklabels=10,100,1000,10000,100000,
width=.9, height=.8, enlarge y limits=0.05, y=0.28cm, bar width=0.25cm,
xlabel=Feature importance (logarithmic scale),

symbolic y coords=area@0.16,area@0.24,area@0.32,area@0.08,area@0.50,value@0.04,blb_sum@0.08,blb_sum@0.50,blb_cnt@0.08,value@0.02,blb_cnt@0.16,area@0.02,blb_cnt@0.50,value@0.08,blb_cnt@0.04,blb_cnt@0.02,value@0.16,blb_sum@0.04,blb_cnt@0.24,blb_sum@0.02,value@total,blb_sum@0.16,area@0.04,blb_sum@0.24,value@0.24,value<EMAIL_ADDRESS>coordinates<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(1388.45,value@total) <EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(182609.1,value@0.50);

GBT top feature importance.

§.§ Discussion

As we can see in Table <ref> and Table <ref>, direct cellularity assessment method slightly outperforms the segmentation-based approach, where the regression module works on top of the segmentation feature extractor. We believe that performance improves due to two main reasons. First, segmentation models were not trained on accurate segmentation masks but rather on approximate masks generated from weakly supervised labels. Second, the cellularity score depends not only on the tumor masks but also on a broader set of features, some of which could be lost during the segmentation step.

While we note the record results of our end-to-end models, we believe that the modular form of the prediction pipeline provides benefits that more than compensate for this small difference in the final score.

The segmentation-based approach has two significant advantages: generalizability and interpretability. In practice, the data used for medical imaging tasks comes from different hospitals and is collected by different hardware. Images may differ in quality, level of noise, color and brightness distributions. In [16], the authors proposed to use segmentation to clean and standardize the data, which helps with overall robustness and performance of various task-specific models.

Better interpretability is achieved by the fact that we can visually verify the quality of the intermediate step, i.e., segmented tumors. Furthermore, the decision trees model allows to estimate the feature importance for every feature based on the information gain. If segmented tumors are correct, and the most informative features make intuitive sense, we obtain additional confidence in our model, which is very important in the medical setting.

§ CONCLUSION

In this paper, we evaluate three automatic methods to assess the cellulalarity of residual breast tumors in H&E stained samples. Our first method leverages the weakly-supervised segmentation masks as inputs for deep CNN. We believe that this method will be more generalizable and robust towards the data acquisition and easier to interpret.

Our second method that leverages feature extraction from the weakly-supervised segmentation mask yields the highest score among the all previously published feature extraction-based methods [3, 29].

Finally, the third method in this study is an end-to-end approach that predicts the cellularity score without any intermediate segmentation step. Although it is attractive and produces the best results it lacks interpretability of the segmentation-based methods and could perform best due to the dataset bias.

The main limitation of this study is the dataset size and the weak labels for the segmentation model. We think that given a bigger dataset and good quality annotations, segmentation-based approach could produce better results that less deviate from the end-to-end trained models.

§ ACKNOWLEDGEMENTS

The work of Sergey Nikolenko was supported by the Russian Foundation for Basic Research grant no. 18-54-74005.
The authors thank the anonymous reviewer #2 for the constructive suggestions that helped to improve the article.

[1]
E. K. V. I. I. A. Buslaev, A. Parinov and A. A. Kalinin.
Albumentations: fast and flexible image augmentations.
ArXiv e-prints, 2018.

[2]
S. Akbar, M. Peikari, S. Salama, S. Nofech-Mozes, and A. L. Martel.
Determining tumor cellularity in digital slides using resnet.
In Medical Imaging 2018: Digital Pathology, volume 10581, page
105810U, 2018.

[3]
S. Akbar, M. Peikari, S. Salama, A. Y. Panah, S. Nofech-Momes, and A. L.
Automated and manual quantification of tumour cellularity in digital
slides for tumour burden assessment.
bioRxiv, page 571190, 2019.

[4]
T. Araújo, G. Aresta, E. Castro, J. Rouco, P. Aguiar, C. Eloy,
A. Polónia, and A. Campilho.
Classification of breast cancer histology images using convolutional
neural networks.
PloS one, 12(6):e0177544, 2017.

[5]
B. E. Bejnordi, M. Veta, P. J. van Diest, B. van Ginneken, N. Karssemeijer,
G. Litjens, J. A. van der Laak, M. Hermsen, Q. F. Manson, M. Balkenhol,
et al.
Diagnostic assessment of deep learning algorithms for detection of
lymph node metastases in women with breast cancer.
Jama, 318(22):2199–2210, 2017.

[6]
M. Berman, A. Rannen Triki, and M. B. Blaschko.
The lovász-softmax loss: A tractable surrogate for the
optimization of the intersection-over-union measure in neural networks.
In The IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), June 2018.

[7]
G. Bradski.
The OpenCV Library.
Dr. Dobb's Journal of Software Tools, 2000.

[8]
Detailed pathology methods for using residual cancer burden.

[9]
T. Ching, D. S. Himmelstein, B. K. Beaulieu-Jones, A. A. Kalinin, B. T. Do,
G. P. Way, E. Ferrero, P.-M. Agapow, M. Zietz, M. M. Hoffman, et al.
Opportunities and obstacles for deep learning in biology and
Journal of The Royal Society Interface, 15(141), 2018.

[10]
L. E. Boucheron, B. Manjunath, and N. Harvey.
Use of imperfectly segmented nuclei in the classification of
histopathology images of breast cancer.
pages 666–669, 03 2010.

[11]
J. G. Elmore, G. M. Longton, P. A. Carney, B. M. Geller, T. Onega, A. N.
Tosteson, H. D. Nelson, M. S. Pepe, K. H. Allison, S. J. Schnitt, et al.
Diagnostic concordance among pathologists interpreting breast biopsy
Jama, 313(11):1122–1132, 2015.

[12]
K. He, X. Zhang, S. Ren, and J. Sun.
Delving deep into rectifiers: Surpassing human-level performance on
imagenet classification.
In Proceedings of the IEEE international conference on computer
vision, pages 1026–1034, 2015.

[13]
K. He, X. Zhang, S. Ren, and J. Sun.
Deep residual learning for image recognition.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 770–778, 2016.

[14]
K. He, X. Zhang, S. Ren, and J. Sun.
Deep residual learning for image recognition.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 770–778, 2016.

[15]
V. Iglovikov and A. Shvets.
Ternausnet: U-net with vgg11 encoder pre-trained on imagenet for
image segmentation.
arXiv preprint arXiv:1801.05746, 2018.

[16]
V. I. Iglovikov, A. Rakhlin, A. A. Kalinin, and A. A. Shvets.
Paediatric bone age assessment using deep convolutional neural
In Deep Learning in Medical Image Analysis and Multimodal
Learning for Clinical Decision Support, pages 300–308. Springer, 2018.

[17]
S. Ioffe and C. Szegedy.
Batch normalization: Accelerating deep network training by reducing
internal covariate shift.
In International Conference on Machine Learning, pages
448–456, 2015.

[18]
P. Izmailov, D. Podoprikhin, T. Garipov, D. Vetrov, and A. G. Wilson.
Averaging weights leads to wider optima and better generalization.
arXiv preprint arXiv:1803.05407, 2018.

[19]
A. A. Kalinin, A. Allyn-Feuer, A. Ade, G.-V. Fon, W. Meixner, D. Dilworth,
S. S. Husain, J. R. de Wett, G. A. Higgins, G. Zheng, et al.
3D shape modeling for cell nuclear morphological analysis and
Scientific Reports, 8, 2018.

[20]
G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu.
Lightgbm: A highly efficient gradient boosting decision tree.
In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus,
S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information
Processing Systems 30, pages 3146–3154. Curran Associates, Inc., 2017.

[21]
D. Komura and S. Ishikawa.
Machine learning methods for histopathological image analysis.
Computational and structural biotechnology journal, 16:34–42,

[22]
Y. LeCun, Y. Bengio, and G. Hinton.
Deep learning.
nature, 521(7553):436, 2015.

[23]
T. Lin, P. Dollár, R. B. Girshick, K. He, B. Hariharan, and S. J.
Feature pyramid networks for object detection.
CoRR, abs/1612.03144, 2016.

[24]
T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, and S. Belongie.
Feature pyramid networks for object detection.
In Proceedings of the IEEE Conference on Computer Vision and
Pattern Recognition, pages 2117–2125, 2017.

[25]
Y. Liu, D. Minh Nguyen, N. Deligiannis, W. Ding, and A. Munteanu.
Hourglass-shapenetwork based semantic segmentation for high
resolution aerial imagery.
Remote Sensing, 9(6):522, 2017.

[26]
K. Lopuhin.
Noaa fisheries steller sea lion population count.
2017, online; accessed April 18, 2019.

[27]
J. S. Meyer, C. Alvarez, C. Milikowski, N. Olson, I. Russo, J. Russo, A. Glass,
B. A. Zehnbauer, K. Lister, and R. Parwaresch.
Breast carcinoma malignancy grading by bloom–richardson system vs
proliferation index: reproducibility of grade and advantages of proliferation
Modern pathology, 18(8):1067, 2005.

[28]
A. Natekin and A. Knoll.
Gradient boosting machines, a tutorial.
Frontiers in neurorobotics, 7:21, 2013.

[29]
M. Peikari, S. Salama, S. Nofech-Mozes, and A. L. Martel.
Automatic cellularity assessment from post-treated breast surgical
Cytometry Part A, 91(11):1078–1087, 2017.

[30]
R. Rajan, A. Poniecka, T. L. Smith, Y. Yang, D. Frye, L. Pusztai, D. J.
Fiterman, E. Gal-Gombos, G. Whitman, R. Rouzier, et al.
Change in tumor cellularity of breast carcinoma after neoadjuvant
chemotherapy as a variable in the pathologic assessment of response.
Cancer: Interdisciplinary International Journal of the American
Cancer Society, 100(7):1365–1373, 2004.

[31]
A. Rakhlin, A. Davydow, and S. Nikolenko.
Land cover classification from satellite imagery with u-net and
lovász-softmax loss.
In The IEEE Conference on Computer Vision and Pattern
Recognition (CVPR) Workshops, June 2018.

[32]
S. Robertson, H. Azizpour, K. Smith, and J. Hartman.
Digital image analysis in breast pathology—from image processing
techniques to artificial intelligence.
Translational Research, 2017.

[33]
O. Ronneberger, P. Fischer, and T. Brox.
U-net: Convolutional networks for biomedical image segmentation.
In International Conference on Medical image computing and
computer-assisted intervention, pages 234–241. Springer, 2015.

[34]
J. Schmidhuber.
Deep learning in neural networks: An overview.
Neural networks, 61:85–117, 2015.

[35]
A. A. Shvets, A. Rakhlin, A. A. Kalinin, and V. I. Iglovikov.
Automatic instrument segmentation in robot-assisted surgery using
deep learning.
In 2018 17th IEEE International Conference on Machine Learning
and Applications (ICMLA), pages 624–628. IEEE, 2018.

[36]
R. L. Siegel, K. D. Miller, and A. Jemal.
Cancer statistics, 2018.
CA: A Cancer Journal for Clinicians, 68(1):7–30, 2018.

[37]
K. Simonyan and A. Zisserman.
Very deep convolutional networks for large-scale image recognition.
arXiv preprint arXiv:1409.1556, 2014.

[38]
F. A. Spanhol, L. S. Oliveira, C. Petitjean, and L. Heutte.
Breast cancer histopathological image classification using
convolutional neural networks.
In Neural Networks (IJCNN), 2016 International Joint Conference
on, pages 2560–2567. IEEE, 2016.

[39]
W. F. Symmans, F. Peintinger, C. Hatzis, R. Rajan, H. Kuerer, V. Valero,
L. Assad, A. Poniecka, B. Hennessy, M. Green, et al.
Measurement of residual breast cancer burden to predict survival
after neoadjuvant chemotherapy.
Journal of Clinical Oncology, 25(28):4414–4422, 2007.

[40]
C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna.
Rethinking the inception architecture for computer vision.
In Proceedings of the IEEE conference on computer vision and
pattern recognition, pages 2818–2826, 2016.

[41]
A. Tiulpin, J. Thevenot, E. Rahtu, P. Lehenkari, and S. Saarakkala.
Automatic knee osteoarthritis diagnosis from plain radiographs: A
deep learning-based approach.
Scientific Reports, 8:1727, 2018.

[42]
J. Tompson, R. Goroshin, A. Jain, Y. LeCun, and C. Bregler.
Efficient object localization using convolutional networks.
In Proceedings of the IEEE Conference on Computer Vision and
Pattern Recognition, pages 648–656, 2015.

[43]
M. Veta, P. J. Van Diest, and J. P. Pluim.
Cutting out the middleman: measuring nuclear area in histopathology
slides without segmentation.
In International Conference on Medical Image Computing and
Computer-Assisted Intervention, pages 632–639. Springer, 2016.

[44]
F. Xing and L. Yang.
Robust nucleus/cell detection and segmentation in digital pathology
and microscopy images: A comprehensive review.
IEEE Reviews in Biomedical Engineering, 9:234–263, 2016.

[45]
R. Zhang, P. Isola, and A. A. Efros.
Colorful image colorization.
Lecture Notes in Computer Science, page 649–666, 2016.