# Grad-CAMO: Learning Interpretable Single-Cell Morphological Profiles from 3D
Cell Painting Images

Vivek Gopalakrishnan
Xellar Biosystems and MIT
<EMAIL_ADDRESS>Jingzhe Ma
Xellar Biosystems
<EMAIL_ADDRESS>Zhiyong Xie
Xellar Biosystems
<EMAIL_ADDRESS>
###### Abstract

Despite their black-box nature, deep learning models are extensively used in
image-based drug discovery to extract feature vectors from single cells in
microscopy images. To better understand how these networks perform
representation learning, we employ visual explainability techniques (_e.g_.
Grad-CAM). Our analyses reveal several mechanisms by which supervised models
cheat, exploiting biologically irrelevant pixels when extracting morphological
features from images, such as noise in the background. This raises doubts
regarding the fidelity of learned single-cell representations and their
relevance when investigating downstream biological questions. To address this
misalignment between researcher expectations and machine behavior, we
introduce Grad-CAMO, a novel single-cell interpretability score for supervised
feature extractors. Grad-CAMO measures the proportion of a model’s attention
that is concentrated on the cell of interest versus the background. This
metric can be assessed per-cell or averaged across a validation set, offering
a tool to audit individual features vectors or guide the improved design of
deep learning architectures. Importantly, Grad-CAMO seamlessly integrates into
existing workflows, requiring no dataset or model modifications, and is
compatible with both 2D and 3D Cell Painting data. Additional results are
available at https://github.com/eigenvivek/Grad-CAMO.

##  1 Introduction

Drug discovery, screening, and development is a lengthy and costly process
with a high failure rate, due in large part to the numerous discrepancies
between cells cultured on 2D surfaces and those inhabiting in vivo cellular
environments [6]. 3D tissue culture techniques, particularly microfluidic-
based organ-on-a-chip platforms [1], enable the development of cellular
disease models that more closely mimic in vivo biology. When combined with
high-content, high-throughput imaging techniques, such as the Cell Painting
assay [2], these advanced tissue culture approaches can be used to measure the
effect of diverse treatments on 3D cellular morphology. However, for cellular
analyses to scale with the immense rate with which modern imaging platforms
acquire data, we also require computational methods that automatically
transform single cells in 3D Cell Painting images into quantitative
representations.

Figure 1: Overview of Grad-CAMO. Given a 3D Cell Painting Z-stack, we first
use Cellpose [12] to segment individual cells. Using segmentation masks, we
create cuboid 3D crops centered on single cells. Extending approaches commonly
used for 2D Cell Painting images (_e.g_. DeepProfiler [9]), we train a 3D
EfficientNet [14] to predict the treatment label of an individual cell from a
crop. During inference, held out cells are passed to the trained network and
activations at intermediate layers are used as single-cell feature vectors in
$\mathbb{R}^{d}$. Visualization of feature vectors using UMAP [7] shows single
cells are highly separable based on learned morphological profiles. However,
interpretability analysis using Grad-CAM demonstrates that deep learning-based
feature extractors do not always pay attention to the cell-of-interest when
forming single-cell morphological profiles (on-target _vs_. off-target). To
quantify the fidelity of learned morphological profiles, we introduce Grad-
CAMO, a single-cell interpretability metric to quantify the level of
confounding in a model’s predictions.

Quantifying the phenotypic effects of experimental perturbations from high-
throughput imaging assays, a feature extraction process known as morphological
profiling, is a necessary and challenging step in the analysis of Cell
Painting data. Whereas traditional methods use handcrafted image processing
algorithms to extract human-designed descriptors of cellular morphology
(_e.g_. CellProfiler [4]), recent methods have taken a deep learning approach,
training a neural network to extract features learned to be relevant from raw
imaging data [9, 8]. While learned morphological profiles have been shown to
enable better performance in downstream analysis tasks, the question of
interpretability limits deep learning-based approaches: how can we ensure that
the morphological profiles extracted by black-box deep learning models capture
biologically relevant information about single cells, rather than simply
exploiting confounding factors present in image data (_e.g_. batch effects) to
minimize training loss?

Previous work has suggested that pretraining supervised feature extractors on
large, diverse sets of 2D Cell Painting images can mitigate the impact of
confounders on learned morphological profiles [9]. Unfortunately, for 3D Cell
Painting, there are no open-source imaging databases as extensive as the JUMP-
Cell Painting or Recursion datasets [5, 13], which contain millions of 2D
microscopy images for pretraining deep learning models. This data limitation
also precludes the development of self-supervised feature extractors (_e.g_.
attention-based mechanisms such as DINO [3]), which have shown promise in 2D
Cell Painting images [10]. As transformer-based models are currently
infeasible in for 3D microscopy data, there remains a need for quantitative
methods to audit the morphological profiles produced by supervised deep
learning models.

To address this need, we propose evaluating learned morphological profiles
with Gradient-weighted Class Activation Mapping (Grad-CAM) [11], a technique
that uses the gradients of a convolutional neural network to localize which
regions of an input image the model paid the most attention to when making its
prediction. When combined with single-cell segmentations masks, we can measure
the proportion of a model’s attention that overlaps with the cell of interest
compared to the background. This interpretability metric, which we term Grad-
CAMO (_i.e_. Grad-CAM Overlap), enables identification of morphological
profiles that faithfully describe biologically relevant components of the
input image, helping to quantify the influence of confounders on the extracted
features.

In our experiments, we demonstrate the utility of Grad-CAMO for 3D Cell
Painting images: using a dataset of single-cell 3D crops extracted from
fluorescence microscopy Z-stacks, we train a 3D convolutional neural network
to predict the treatment label of each cell. Intermediate activations are
extracted and used as single-cell features. Using Grad-CAMO, we find that only
30% of learned morphological profiles have Grad-CAM localization maps that
meaningfully overlap with the cell’s segmentation map. Visualization of the
localization maps demonstrates that supervised feature extractors can cheat by
exploiting irrelevant, non-biological information in microscopy data.

#### Contributions.

We introduce Grad-CAMO, an easily computed interpretability metric that
quantifies to what extent a single-cell morphological profile represents the
cell of interest. Like the Grad-CAM localization map, Grad-CAMO is data- and
model-agnostic: it can be computed for any convolutional feature extractor
without modification to the imaging workflow or deep learning architecture.
However, unlike previous interpretability techniques, which only provide a
qualitative assessment of fidelity, Grad-CAMO provides a quantitative measure
of the biological relevance of the morphological profile for each single cell.
This highlights the utility of our metric for screening computational
pipelines that process millions of cells from high-throughput imaging data.
Although this work primarily focused on 3D images, Grad-CAMO is also
extensible to 2D Cell Painting images and can be used to explain descriptions
made by existing models employed in academia and industry. The integration of
Grad-CAMO into a traditional single-cell feature extraction pipeline is shown
in Figure 1.

##  2 Preliminaries

Let $\mathbf{I}:\mathbb{R}^{3}\to\mathbb{R}^{C}$ represent a Z-stack
comprising thousands of cells, where $\mathbf{I}$ is a function mapping 3D
coordinates to pixel intensities in a 3D Cell Painting image with $C$
fluorescent channels. Let $\mathbf{V}_{i}\subseteq\mathbf{I}$ be a single-cell
3D crop of the Z-stack centered on the $i$-th cell. Finally, let
$Y_{i}\in\mathcal{Y}=\\{1,\dots,K\\}$ be a label for the $i$-th cell,
corresponding to a categorical dimension (_e.g_. the treatment dosage) along
which we wish to extract morphological profiles for some downstream task. We
assume there is at least one Z-stack for each categorical label and that
images come from wells wherein each cell received the same treatment.

###  2.1 Supervised Feature Extraction

To extract morphological profiles using supervised learning, we first train a
convolutional neural network $f_{\theta}$ to predict $Y_{i}$ from
$\mathbf{V}_{i}$ by minimizing a classification loss with respect to weights
$\theta$ over some subset of cells in the Z-stacks. Then, at inference time,
cells not seen by the network during training are processed by the network and
activations at an intermediate layer are extracted as learned morphological
profiles. To use explicit notation, let us decouple $f_{\theta}$ into a
convolutional backbone $g_{\varphi}$ and a classification head $h_{w}$ such
that

$\hat{Y}_{i}=f_{\theta}(\mathbf{V}_{i})=(h_{w}\circ
g_{\varphi})(\mathbf{V}_{i})\,,$ (1)

where $\hat{Y}_{i}$ is the predicted treatment label for cell
$\mathbf{V}_{i}$. Then, for a previously unseen cell $\mathbf{V}$, the
intermediate activation

$\mathbf{A}=g_{\varphi}(\mathbf{V})\in\mathbb{R}^{C^{\prime}\times H\times
W\times D}$ (2)

is a representation of the single-cell crop that the network has learned is
useful for predicting the treatment label. This representation can then be
pooled and reshaped into a feature vector in $\mathbb{R}^{d}$. Representing 3D
Cell Painting images as a set of single-cell morphological profiles enables
numerous downstream machine learning tasks such as clustering and
classification.

###  2.2 Grad-CAM

While a neural network given sufficient data will eventually learn to
accurately estimate the treatment label of a single-cell crop, it is not
obvious how the network arrived at its prediction. This calls into the
question the fidelity of the intermediate activation $\mathbf{A}$, as it could
ostensibly arise from a spurious correlation or confounding variable latent in
the input data. To visually explain the decisions made by a convolutional
neural network, we can use Gradient-weighted Class Activation Mapping (Grad-
CAM) [11], an interpretability technique that exploits the model’s gradients
to produce a localization map highlighting which regions in the input image
led the model to predict a particular label. Specifically, for a cell that the
model predicts received treatment $k\in\mathcal{Y}$, the Grad-CAM localization
map is a linear combination of the channels in $\mathbf{A}$, calculated as

$\tilde{\mathbf{G}}=\mathrm{ReLU}\left(\frac{1}{C^{\prime}}\sum_{c=1}^{C^{\prime}}w_{c}^{(k)}\mathbf{A}_{c}\right)\in\mathbb{R}^{H\times
W\times D}\,,$ (3)

where the weight $w_{c}^{(k)}$ is computed with per-channel average pooling:

$w_{c}^{(k)}=\frac{1}{HWD}\sum_{x=1}^{H}\sum_{y=1}^{W}\sum_{z=1}^{D}\frac{\partial
Y^{(k)}}{\partial\mathbf{A}_{c}(x,y,z)}\,,$ (4)

where $Y^{(k)}$ is the model’s score for treatment label $k$ (_i.e_. the
output prior to applying the softmax). This weight captures the importance of
the $c$-th channel to the network’s prediction, and the ReLU in Eq. 3 serves
to preserve only those regions that positively influence this prediction.

Note that $\tilde{\mathbf{G}}$ only provides a coarse localization since it is
the same shape as $\mathbf{A}$, and intermediate activations at bottleneck
layers in convolutional neural networks typically have a much smaller spatial
dimensions than the input. Therefore, we upsample $\tilde{\mathbf{G}}$ to
match the original dimension of $\mathbf{V}$ using trilinear interpolation and
refer to the rescaled Grad-CAM localization map as
$\mathbf{G}=\mathrm{Trilinear}(\tilde{\mathbf{G}})$.

##  3 Methods

###  3.1 Grad-CAMO

While Grad-CAM can be used to visually inspect the outputs of deep networks,
it is infeasible to manually inspect every localization map for the thousands
of cells in a Z-stack. Therefore, we introduce Grad-CAMO (short for Grad-CAM
Overlap), a single number score to quantify the concentration of the
localization map within the cell of interest. Let $\mathbf{M}$ be a binary
segmentation mask with the same shape as $\mathbf{V}$ that denotes the
occupancy of the central cell within the volume. Let $\mathbf{G}$ be the
localization map produced upsampled to match the dimensions of $\mathbf{V}$.
Then, we can compute the Grad-CAMO score as

$\mathrm{Grad\text{-}CAMO}(\mathbf{G},\mathbf{M})=\frac{\|\mathrm{vec}(\mathbf{G}\odot\mathbf{M})\|_{1}}{\
|\mathrm{vec}(\mathbf{G})\|_{1}}\,,$ (5)

where $\odot$ represents element-wise matrix multiplication.

Grad-CAMO is a model- and data-agnostic score, able to be automatically
computed for any convolutional neural network. Additionally, although not
explored in this work, Grad-CAMO is extensible to vision transformer-based
feature extractors. This can be accomplished simply by replacing the Grad-CAM
localization map $\mathbf{G}$ with the self-attention map produced by vision
transformer models.

Figure 2: Automatic segmentations produced by Cellpose. Zoomed-in image
patches are shown to demonstrate the accuracy of segmentations produced by
Cellpose with minimal hyperparameter tuning. 2D segmentation masks were
stitched across the $z$-axis of the 3D Cell Painting images to form 3D
segmentation masks.

###  3.2 Model Auditing

The Grad-CAMO score ranges from $[0,1]$, with higher scores denoting
morphological profiles that more closely relate to the cell of interest.
Because this score can be automatically evaluated for every cell independent
of model architecture, Grad-CAMO can be used to audit supervised feature
extractors, either by filtering cells with biologically irrelevant
morphological profiles or by providing a heuristic to guide improvements to
the design of the feature extractor. For example, let
$\mathcal{D}_{\mathrm{val}}=\\{\mathbf{V}_{1},\dots,\mathbf{V}_{N}\\}$ be a
validation set of $N$ held-out single-cell crops. For each cell, we can
compute a Grad-CAMO score $\\{s_{1},\dots,s_{N}\\}$. Optimal design choices
(_e.g_. what layer to extract features from, what pooling function to use,
_etc_.) can be empirically determined by maximizing a central measure of Grad-
CAMO over the validation dataset, _e.g_.
$\hat{s}=\frac{1}{N}\sum_{i=1}^{N}s_{i}$. Thus, Grad-CAMO a useful metric for
hyperparameter tuning.

###  3.3 Implementation Details

#### Single-cell segmentation.

We segment individual cells in 3D Cell Painting images using the cyto model in
Cellpose [12], which only uses the nucleus and cytoplasm channels in the
Z-stack. We do not use the built-in 3D segmentation model, but rather achieve
higher performance by segmenting every slice in the Z-stack with the 2D model
and combining all slices into a 3D segmentation mask using a stitch threshold
of 0.05. Additional non-default hyperparameters used include a flow threshold
of 0.75, an anisotropy factor of $1.71$ (_i.e_. the ratio of between-plane
pixel spacing to in-plane pixel spacing), and a minimum pixel size of 30
pixels per 2D segmentation mask. A 2D slice of the 3D segmentation masks
produced by Cellpose is shown in Figure 2.

Using 3D segmentation masks, $128\times 128\times Z$ single-cell crops were
extracted from the Z-stack images, where $Z=21$ is the total number of slices
in the Z-stack. This was done as nearly all cells spanned the entire $z$-axis.
For each crop, the cell center was determined as the average of the minimum
and maximum pixel index in the $x$\- and $y$-axes of the cell’s 3D
segmentation mask. For preprocessing, channels in each single-cell crop were
first rescaled within $[0,1]$, then normalized by the mean and variance of
rescaled pixel intensities in the training set. A 2D slice of a single-cell
crop produced by our segmentation and preprocessing pipeline is shown in
Figure 3.

Figure 3: Examples of multi-channel single-cell crops at different dosages.
For every cell segmented by Cellpose, a $128\times 128\times 21$ crop was
extracted from the Z-stack and preprocessed to standardize pixel intensity
ranges across samples.

#### Supervised learning.

We trained a 3D EfficientNet-B0 to predict the treatment label of single-cell
crops.111We used an open-source implementation available at
https://github.com/shijianjian/EfficientNet-PyTorch-3D. To ensure that single-
cell crops were divisible by the number of pooling layers in the model, we
upsampled all 3D volumes from their original size of $128\times 128\times 21$
voxels to $224\times 224\times 64$ voxels using trilinear interpolation. After
normalization, all single-cell crops were augmented as follows:

  1. 1.

Random $x$\- and $y$-axis flips: in our implementation, the $z$-axis is never
reversed.

  2. 2.

Random brightness adjustment: add a random constant intensity from
$\mathcal{U}[0.5,1.25]$ to all pixel intensities.

  3. 3.

Random gamma adjustment: raise all pixel intensities to a random power sampled
from $\mathcal{U}[0.5,1.5]$.

We do not include random crop and resize augmentations as they alter cell
size, a biologically relevant feature that can provide information about the
health of a cell under different treatments. All augmentations were applied
randomly with probability $p=0.5$. During training, we minimized cross entropy
loss using the Adam optimizer with a learning rate of $3\times 10^{-4}$. All
models were trained on a single NVIDIA A6000 GPU with a batch size of 8 single
cells.

#### Correcting for batch effects.

We employ the whitening transform fit to the control cells as our batch effect
correction strategy. In statistical terms, if $\mathbf{x}\in\mathbb{R}^{d}$ is
a random feature vector with mean $\mathbf{0}$ and covariance
$\mathbf{\Sigma}$, the matrix $\mathbf{W}\in\mathbb{R}^{d\times d}$ performs
the whitening transform if $\mathbf{W}\mathbf{x}$ has a covariance matrix
equal to identity. That is, the whitening transform maps samples from a
particular distribution to isotropic Gaussian noise. Given a (centered) data
matrix $\mathbf{X}\in\mathbb{R}^{N\times d}$ whose rows are feature vectors
from $N$ control cells, we estimate
$\mathbf{W}=\hat{\mathbf{\Sigma}}^{-\frac{1}{2}}$ via an eigendecomposition of
the sample covariance matrix
$\hat{\mathbf{\Sigma}}=\frac{1}{N}\mathbf{X}^{T}\mathbf{X}$ [9]. Cells in the
control group are used to estimate the whitening transform as they should have
identical phenotypic profiles, and thus, any differences in their features are
sampled from the underlying noise distribution. Given feature vectors
$\mathbf{Y}\in\mathbb{R}^{M\times d}$ from other cells, we use the batch-
corrected features $\mathbf{W}\mathbf{Y}^{T}$.

##  4 Results

###  4.1 Dataset

Our 3D Cell Painting dataset was derived from a drug-induced liver injury
study of Hep G2 cells treated with Stavudine at six different concentrations.
Cells were seeded in a 3D growth medium and cultured for 14 days. Two wells
were cultured at each concentration, and four distinct sites were imaged per
well. Cells were stained with six fluorescent dyes including Hoescht
(nucleus), Concanavalin (endoplasmic reticulum), SYTO14 (nucleoli and
cytoplasmic RNA), Wheat-Germ Agglutinin (cell membrane and Golgi apparatus),
Phalloidin (F-actin cytoskeleton), and MitoTracker (mitochondria) [2]. 3D Cell
Painting data were acquired using six channels at 20$\times$ magnification and
an anisotropic voxel spacing of $$0.3505\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$0.3505\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$0.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$$. For each volume, 21
slices were acquired along the $z$-axis.

Table 1: 3D Cell Painting dataset. Six treatment groups were cultured with two biological replicates and four technical replicates. Cell counts were estimated using segmentations from Cellpose. Treatment Dose | Well | Number of Cells
---|---|---
Control | D02, D03 | 9,313
0.06$\times$ | B02, B03 | 4,119
0.32$\times$ | C02, C03 | 10,182
1.60$\times$ | G02, G03 | 9,468
8.00$\times$ | F02, F03 | 9,813
40.0$\times$ | E02, E03 | 9,416
Figure 4: (A) UMAP and Grad-CAM. UMAP embedding of features extracted from the
single cells in the test set shows that the different treatment labels are
highly separable: ( ) Control, ( ) $0.06\times$, ( ) $0.32\times$, ( )
$1.60\times$, ( ) $8.00\times$, and ( ) $40.0\times$. Correction for batch
effects is accomplished via a whitening transform. However, cells also cluster
by site ($\mathbin{\vbox{\hbox{\scalebox{1.5}{$\bullet$}}}}$,
$\mathbf{\times}$, $\mathbin{\vbox{\hbox{\scalebox{1.0}{$\blacksquare$}}}}$,
$\mathbf{+}$), demonstrating the effect of confounding on learned
morphological profiles. (B) Where is the model looking? Using Grad-CAM, we
identify three patterns in model attention during deep morphological
profiling: concentrating on the central cell, a neighboring cell, or the
background. In these localization maps, red denotes higher attention. For
visualization purposes, only render the central slice of the Z-stack and
overlay associated slice in the Grad-CAM localization map.

Segmentation with Cellpose enumerated 52,311 cells across all biological and
technical replicates (Table 1), which were subsequently processed into a
dataset of single-cell crops. Partial cells on the boundary of a Z-stack were
removed from the dataset and are not included in the total cell counts. To
mitigate the potential for data leakage between the training and testing
datasets, we used a leave-wells-out validation strategy, training all models
on cells from wells in the second column of the plate and evaluated on cells
in the third column. A subset of cells within the training set were used as a
held-out validation set for an early stopping criterion.

Figure 5: (A) Example Grad-CAMO scores. Grad-CAMO is calculated as the
proportion of the model’s Grad-CAM localization map that lies within the
segmentation mask of the central cell. (B) Grad-CAMO scores computed over the
entire testing set. Distributions are grouped by treatment dose and which site
in the well was imaged (_i.e_. technical replicate): ( ) Site 1, ( ) Site 2, (
) Site 3, ( ) Site 4.

###  4.2 Visualizing Learned Morphological Profiles

A 3D EfficientNet-B0 was trained to predict the treatment label of single-cell
crops for ten epochs (or roughly two hours). This short training time was
sufficient for both training and validation accuracy to plateau at 92% and
89%, respectively, demonstrating the simplicity of this classification task
for convolutional models. At the end of each epoch, morphological profiles
were extracted from single cells in the test set and visualized using UMAP
(Figure 4A). Prior to dimensionality reduction, a whitening transform centered
on cells in the control group was performed on extracted features in an
attempt to eliminate batch effects. Visualization of the learned morphological
profiles following batch correction demonstrates that cells in the test set
are also highly separable. That is, deep learning strategies extrapolated from
2D use-cases are also able to learn representations of 3D Cell Painting images
that are discriminatory based on a categorical variable (here, the model
achieves 85% accuracy on the test set). However, high classification accuracy
belies a strong source of confounding in the extracted features: cells in the
UMAP cluster not only by the treatment dose, but also by which of the four
sites the cell was imaged from (_i.e_. cells also cluster by technical
replicate). While the whitening transform removes this source of variation for
control cells by mapping their morphological profiles to a standard Gaussian
distribution, this batch correction method is insufficient to recover faithful
morphological representations for the other treatment doses.

Visualizing Grad-CAM localization maps reveals three attention patterns when
the model is predicting the treatment of a previously unseen cell: (1) the
model concentrates on the central cell of interest, as intended; (2) the model
concentrates on the morphology of a neighboring cell; and (3) the model
concentrates on some biologically irrelevant noise in the background of the
image. These three modes are visualized in Figure 4B.

This experiment illustrates two mechanisms by which the morphological profiles
estimated by supervised machine learning models are misaligned with the
expectations of human practitioners. First, we see that models achieve high
accuracy by sometimes classifying neighboring cells, perhaps because the
neighbors in these cases closely resemble a cell in the training set. However,
this is incongruous with our assumption that single-cell feature vectors are
describing the central cell in a 3D crop. Second, this visualization also
demonstrates that supervised feature extractors can simply cheat by exploiting
non-biological information in microscopy data. In the examples in Figure 4B,
there are no human-discernible features in the background of the image that
would enable the practitioner to determine the treatment received by these
cells. While the UMAP embedding demonstrates highly distinct clusters, visual
interpretation of these morphological profiles show that deep learning models
do not always capture the biological information we expect them to.

###  4.3 Computing Grad-CAMO Scores

Visualizing Grad-CAM localization maps helps diagnose failure modes for
supervised feature extractors. However, manually inspecting every
morphological profile in this manner would be too labor-intensive. To
efficiently evaluate our trained model, we instead compute the Grad-CAMO score
for every cell in the test set.

An illustration of how Grad-CAMO is calculated is shown in Figure 5A. In the
first example, the Grad-CAM localization map is centered on the cell of
interest as desired, resulting in a Grad-CAMO score of 0.94. In the second
example, the model concentrates on the unusual cellular morphology, noting the
absence of a round nucleus. However, the model also concentrates on
neighboring cells, resulting in a Grad-CAMO score of 0.24. In the final
example, the model only concentrates on a neighboring cell and receives a
Grad-CAMO score of 0.11. Cells in the test set had an average Grad-CAMO score
of $0.26\pm 0.34$. Plots of the distribution of Grad-CAMO scores, grouped by
treatment dose and site, show that batch effects present in certain treatment
dosages confound the morphological profiles extracted by this model (Figure
5B).

While the model achieves high classification accuracy on the test set and
feature vectors visualized via UMAP cluster neatly, analysis with Grad-CAMO
demonstrates that these single-cell morphological profiles have low
interpretability and poorly represent the central cell in each crop. In fact,
we find that only 30% of learned morphological profiles have Grad-CAM
localization maps that significantly overlap with the cell’s segmentation map
(using a Grad-CAMO score of $0.25$ as a cutoff).

Filtering morphological profiles with low Grad-CAMO scores is a potential
strategy for removing additional batch effects. For example, see the UMAP in
Figure 1, which only plots morphological profiles with a Grad-CAMO score
greater than $0.25$. However, this removes roughly 13,000 cells from the
dataset, reducing the power of the study. Alternatively, Grad-CAMO could be
used to evaluate different design choices for a feature extraction pipeline.
After determination of an optimal model architecture, low Grad-CAMO can be
used to filter outlier cells with unrepresentative morphological profiles.
From a performance point of view, Grad-CAMO is a fast and memory-efficient
operation that can be calculated in parallel with the feature extraction step,
increasing computational overhead by from $22\text{\,}\mathrm{min}$ to
$34\text{\,}\mathrm{min}$ for 23,000 cells on a single GPU. Computing Grad-
CAMO on an entire dataset can also trivially be parallelized across multiple
GPUs.

##  5 Discussion

Given the critical importance of single-cell feature extraction for numerous
downstream biological tasks (_e.g_. measuring dose response, elucidating
mechanisms of action, or predicting clinical outcomes, _etc_.), ensuring the
fidelity and interpretability of morphological profiles is paramount. Deep
learning-based feature extraction has many advantages over classical image
processing pipelines using manually engineered features. For example, deep
models are more computationally efficient, as they execute on the GPU, and are
therefore better suited for high-throughput imaging data. Most importantly,
they can learn rich feature representations that describe phenotypes not
captured by hand-crafted metrics. However, the increased expressiveness of
black-box models comes at a cost: our experiments identify multiple mechanisms
by which the features extracted by neural networks are misaligned with the
expectations of human analysts, including describing the morphologies of
neighboring cells or simply cheating by exploiting imperceptible noise
patterns in the image background.

Grad-CAMO, our proposed explainability metric, can distill the degree of
biological relevance of a morphological profile using a simple and
interpretable score. This score can be evaluated per-sample, to filter
unrepresentative morphological profiles, or across an entire dataset, to
measure the impact of various architectural design choices on feature
extraction pipelines. Maximizing Grad-CAMO serves to optimize the fidelity of
morphological profiles to cells of interest while maintaining the superior
expressiveness of deep learning models. To this end, Grad-CAMO can be used to
align the capabilities of powerful supervised feature extractors with our
expectations of how they perform representation learning.

By improving the explainability and interpretability of supervised deep
learning models, Grad-CAMO has the potential to improve the utility of single-
cell morphological profiles. To conclude, we outline future directions for
Grad-CAMO that we are curious to explore.

#### Grad-CAMO as a heuristic for hyperparameter tuning.

The design space of feature extraction pipelines is incredibly large,
comprising multiple data preprocessing strategies, augmentation options,
convolutional backbones, layers for feature extraction, choices of pooling
layers, _etc_. Oftentimes, design choices are made arbitrarily, and in the
absence of a specified downstream task for which there are ground truth
labels, evaluation of these choices is difficult. To overcome this challenge,
Grad-CAMO can serve as a heuristic to evaluate the impact of different
architectural choices.

#### Grad-CAMO as a regularizer during training.

It is desirable that the morphological profile extracted by a deep learning
model represent the central cell in a 3D crop, instead of a neighboring cell
or the background. This can be enforced during training by using Grad-CAMO as
a regularizer. That is, rather than throwing away single-cell segmentation
masks after data preprocessing, these labels can be used to compute Grad-CAMO
during training. In addition to a classification loss, we can also minimize
negative Grad-CAMO computed across all cells. This form of regularization will
condition intermediate activations of neural networks to focus on the central
cell.

#### Cell Painting inter-channel correlation.

A particularly interesting of Grad-CAMO left unexplored in this work is its
extension to vision transformers-based models. Unlike convolutional
architectures, vision transformers have a built-in mechanism that enables
visual explainability: the self-attention map. The self-attention map can be
used in place of the Grad-CAM localization map to calculate Grad-CAMO for
vision transformers. The current lack of open-source datasets of 3D Cell
Painting images hinders the development of fully self-supervised 3D vision
transformers for morphological profiling. However, supervised transformer-
based models can also exploit the insights gained by through Grad-CAMO to
quantitative measure the quality of morphological profiles obtained by this
increasingly important class models.

## References

  * Bai et al. [2024] Haiqing Bai, Kristen N Peters Olson, Ming Pan, Thomas Marshall, Hardeep Singh, Jingzhe Ma, Paige Gilbride, Yu-Chieh Yuan, Jenna McCormack, Longlong Si, et al.  Rapid prototyping of thermoplastic microfluidic 3D cell culture devices by creating regional hydrophilicity discrepancy.  _Advanced Science_ , 11(7):2304332, 2024.
  * Bray et al. [2016] Mark-Anthony Bray, Shantanu Singh, Han Han, Chadwick T Davis, Blake Borgeson, Cathy Hartland, Maria Kost-Alimova, Sigrun M Gustafsdottir, Christopher C Gibson, and Anne E Carpenter.  Cell Painting, a high-content image-based assay for morphological profiling using multiplexed fluorescent dyes.  _Nature Protocols_ , 11(9):1757–1774, 2016.
  * Caron et al. [2021] Mathilde Caron, Hugo Touvron, Ishan Misra, Hervé Jégou, Julien Mairal, Piotr Bojanowski, and Armand Joulin.  Emerging properties in self-supervised vision transformers.  In _Proceedings of the IEEE International Conference on Computer Vision_ , 2021.
  * Carpenter et al. [2006] Anne E Carpenter, Thouis R Jones, Michael R Lamprecht, Colin Clarke, In Han Kang, Ola Friman, David A Guertin, Joo Han Chang, Robert A Lindquist, Jason Moffat, et al.  CellProfiler: image analysis software for identifying and quantifying cell phenotypes.  _Genome Biology_ , 7:1–11, 2006.
  * Chandrasekaran et al. [2023] SN Chandrasekaran et al.  JUMP Cell Painting dataset: morphological impact of 136,000 chemical and genetic perturbations en. pages: 2023.03. 23.534023 section: New results. mar. 2023, 2023.
  * Langhans [2018] Sigrid A Langhans.  Three-dimensional in vitro cell culture models in drug discovery and drug repositioning.  _Frontiers in Pharmacology_ , 9:334617, 2018.
  * McInnes et al. [2018] Leland McInnes, John Healy, and James Melville.  UMAP: Uniform manifold approximation and projection for dimension reduction.  _arXiv preprint arXiv:1802.03426_ , 2018.
  * Moen et al. [2019] Erick Moen, Dylan Bannon, Takamasa Kudo, William Graf, Markus Covert, and David Van Valen.  Deep learning for cellular image analysis.  _Nature Methods_ , 16(12):1233–1246, 2019.
  * Moshkov et al. [2024] Nikita Moshkov, Michael Bornholdt, Santiago Benoit, Matthew Smith, Claire McQuin, Allen Goodman, Rebecca A Senft, Yu Han, Mehrtash Babadi, Peter Horvath, et al.  Learning representations for image-based profiling of perturbations.  _Nature Communications_ , 15(1):1594, 2024.
  * Pfaendler et al. [2023] Ramon Pfaendler, Jacob Hanimann, Sohyon Lee, and Berend Snijder.  Self-supervised vision transformers accurately decode cellular state heterogeneity.  _bioRxiv_ , pages 2023–01, 2023.
  * Selvaraju et al. [2017] Ramprasaath R Selvaraju, Michael Cogswell, Abhishek Das, Ramakrishna Vedantam, Devi Parikh, and Dhruv Batra.  Grad-CAM: Visual explanations from deep networks via gradient-based localization.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 618–626, 2017.
  * Stringer et al. [2021] Carsen Stringer, Tim Wang, Michalis Michaelos, and Marius Pachitariu.  Cellpose: a generalist algorithm for cellular segmentation.  _Nature Methods_ , 18(1):100–106, 2021.
  * Sypetkowski et al. [2023] Maciej Sypetkowski, Morteza Rezanejad, Saber Saberian, Oren Kraus, John Urbanik, James Taylor, Ben Mabey, Mason Victors, Jason Yosinski, Alborz Rezazadeh Sereshkeh, et al.  RxRx1: A dataset for evaluating experimental batch correction methods.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 4284–4293, 2023.
  * Tan and Le [2019] Mingxing Tan and Quoc Le.  EfficientNet: Rethinking model scaling for convolutional neural networks.  In _International Conference on Machine Learning_ , pages 6105–6114. PMLR, 2019.


$\hat{\mathcal{O}^{\prime}}\Subset\hat{\mathcal{O}}_{0}^{\prime}$ open and
bounded. Moreover by $|\partial\mathcal{E}|_{\mathbb{G}}$ is a Radon measure
then by standard approximation argument we can rewrite $\eqref{cit109}$ with
$\hat{\mathcal{O}^{\prime}}$ instead of $\hat{\mathcal{O}}_{0}^{\prime}$. In
particular using $\eqref{cit106}$ we have that a.e. $c\in(0,1)$

$\liminf_{\alpha\to
0^{+}}|\partial\mathcal{E}_{\alpha,c}|_{\mathbb{G}}(\mathcal{A})=|\partial\mathcal{E}|_{\mathbb{G}}(\mathcal{A}).$

Hence there exists $\hat{c}\in(0,1)$ and a positive sequence
$(\alpha_{h})_{h}$ converging to $0$ such that $\eqref{cit101}$ follows.

Finally because
$\|\phi\|_{\mathcal{L}^{\infty}(\mathcal{O}^{\prime},\mathbb{R})}=M$ we get

$\partial\mathcal{E}\cap\mathcal{A}_{1}=\partial\mathcal{E}\cap\mathcal{A}_{2}=\emptyset$

where
$\mathcal{A}_{1}:=\\{p\in\hat{\mathcal{O}^{\prime}}\cdot\mathbb{V}\,:\,p_{1}\leq-2M\\}$
and
$\mathcal{A}_{2}:=\\{p\in\hat{\mathcal{O}^{\prime}}\cdot\mathbb{V}\,:\,p_{1}\geq
2M\\}$.

Hence since Proposition 5.5 and well-know $\mathbb{G}$-perimeter properties it
follows that

$\displaystyle\frac{1}{C(\mathbb{W},\mathbb{V})}\int_{\mathcal{O}^{\prime}}\sqrt{1+|\nabla^{\hat{\phi}}\hat{\phi}|_{\mathbb{R}^{m-1}}^{2}}\,d\mathcal{L}^{m+n-1}$
$\displaystyle=|\partial\mathcal{E}|_{\mathbb{G}}(\hat{\mathcal{O}^{\prime}}\cdot\mathbb{V})$
$\displaystyle=|\partial\mathcal{E}|_{\mathbb{G}}(\mathcal{A}_{1})+|\partial\mathcal{E}|_{\mathbb{G}}(\mathcal{A})+|\partial\mathcal{E}|(\mathcal{A}_{2})$
$\displaystyle=|\partial\mathcal{E}|_{\mathbb{G}}(\partial\mathcal{E}\cap\mathcal{A}_{1})+|\partial\mathcal{E}|_{\mathbb{G}}(\mathcal{A})+|\partial\mathcal{E}|_{\mathbb{G}}(\partial\mathcal{E}\cap\mathcal{A}_{2})$
$\displaystyle=|\partial\mathcal{E}|_{\mathbb{G}}(\mathcal{A}).$

In the same way, using also $\eqref{cit76}$, $\eqref{cit79}$ and
$\eqref{90cit}$

$|\partial\mathcal{E}_{\alpha_{h},c}|_{\mathbb{G}}(\mathcal{A})=|\partial\mathcal{E}_{\hat{\phi}_{h}}|_{\mathbb{G}}(\mathcal{A})$

and

$|\partial\mathcal{E}_{\hat{\phi}_{h}}|_{\mathbb{G}}(\hat{\mathcal{O}^{\prime}}\cdot\mathbb{V})=\frac{1}{C(\mathbb{W},\mathbb{V})}\int_{\mathcal{O}^{\prime}}\sqrt{1+|\nabla^{\hat{\phi}_{h}}\hat{\phi}_{h}|_{\mathbb{R}^{m-1}}^{2}}\,\,d\mathcal{L}^{m+n-1}$

with $\hat{\phi}_{h}=\hat{\phi}_{\alpha_{h}}$. Now by the last three
equalities and $\eqref{cit101}$, $\eqref{citt98}$ holds. Consequently
$\eqref{citt97}$ is true and the proof of theorem is complete. ∎

##  7\. Intrinsic Lipschitz Function and Distributional Solution

The main results in this section are Theorem 7.1 and Theorem 7.2. We show that
a map $\hat{\phi}=(\phi,0\dots,0):\hat{\mathcal{O}}\to\mathbb{V}$ is a locally
intrinsic Lipschitz function if and only if $\phi$ is locally $1/2$-Hölder
continuous and it is also distributional solution of system

$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w\quad\mbox{in
}\mathcal{O}.$

More precisely, we generalize Theorem 1.1 on [5] valid in Heisenberg groups.

######  Theorem 7.1.

Let $\mathbb{G}:=(\mathbb{R}^{m+n},\cdot,\delta_{\lambda})$ be a Carnot group
of step $2$ and $\mathbb{V}$, $\mathbb{W}$ the complementary subgroups defined
in (8). Let $\hat{\phi}:\hat{\mathcal{O}}\to\mathbb{V}$ be a locally intrinsic
Lipschitz function, where $\hat{\mathcal{O}}$ is an open subset of
$\mathbb{W}$ and $\phi:\mathcal{O}\to\mathbb{R}$ is the map associated to
$\hat{\phi}$ as in $\eqref{phipsi}$.

Then if
$(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi)\in\mathcal{L}_{loc}^{\infty}(\mathcal{O},\mathbb{R}^{m-1})$,
we have that $\phi$ is a distributional solution of
$(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi)=w$ for a suitable
$w\in\mathcal{L}^{\infty}_{loc}(\mathcal{O},\mathbb{R}^{m-1})$ such that
$w(a)=(D^{\phi}_{2}\phi(a),\dots,D^{\phi}_{m}\phi(a))$
$\mathcal{L}^{m+n-1}$-a.e. in $\mathcal{O}$.

###### Proof.

By Theorem 6.1 we know that there exist
$(\phi_{h})_{h\in\mathbb{N}}\subset\mathbb{C}^{\infty}(\mathcal{O},\mathbb{R})$
such that for all $\mathcal{O}^{\prime}\Subset\mathcal{O}$ we have

  1. (1)

$\phi_{h}$ uniformly converges to $\phi$ in $\mathcal{O}^{\prime}$

  2. (2)

$|(D^{\phi_{h}}_{2}\phi_{h}(a),\dots,D^{\phi_{h}}_{m}\phi_{h}(a))|_{\mathbb{R}^{m-1}}\leq\|(D^{\phi}_{2}\phi(a),\dots,D^{\phi}_{m}\phi(a))\|_{\mathcal{L}^{\infty}(\mathcal{O}^{\prime},\mathbb{R}^{m-1})}$
$\forall a\in\mathcal{O}^{\prime}$, $h\in\mathbb{N}$

  3. (3)

for each $j=2,\dots,m$, we have $D^{\phi_{h}}_{j}\phi_{h}(a)\to
D^{\phi}_{j}\phi(a)$ a.e. $a\in\mathcal{O}^{\prime}$

Let
$w_{h}=(w_{2,h},\dots,w_{m,h}):=(D_{2}^{\phi_{h}}\phi_{h},\dots,D_{m}^{\phi_{h}}\phi_{h})$.
Observe that for each $h\in\mathbb{N}$ and
$\zeta\in\mathbb{C}^{\infty}_{c}(\mathcal{O},\mathbb{R})$

$\int_{\mathcal{O}}\phi_{h}\left(\,X_{j}\zeta+\phi_{h}\sum_{s=1}^{n}b_{j1}^{(s)}Y_{s}\zeta\right)\,d\mathcal{L}^{m+n-1}=-\int_{\mathcal{O}}w_{j,h}\zeta\,d\mathcal{L}^{m+n-1}$

for all $j=2,\dots,m$. Getting to the limit for $h\to\infty$ we have

$\int_{\mathcal{O}}\phi\left(\,X_{j}\zeta+\phi\sum_{s=1}^{n}b_{j1}^{(s)}Y_{s}\zeta\right)\,d\mathcal{L}^{m+n-1}=-\int_{\mathcal{O}}w_{j}\zeta\,d\mathcal{L}^{m+n-1}.$

Hence $\phi$ is a distributional solution of
$(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi)=w$ in $\mathcal{O}$. ∎

######  Theorem 7.2.

Let $\mathbb{G}:=(\mathbb{R}^{m+n},\cdot,\delta_{\lambda})$ be a Carnot group
of step $2$ and $\mathbb{V}$, $\mathbb{W}$ the complementary subgroups defined
in (8). Let $\hat{\phi}:\hat{\mathcal{O}}\to\mathbb{V}$ be a continuous map
where $\hat{\mathcal{O}}$ is open in $\mathbb{W}$ and
$\phi:\mathcal{O}\to\mathbb{R}$ is the map associated to $\hat{\phi}$ as in
$\eqref{phipsi}$. We also assume that

  1. (1)

$\phi$ is a continuous distributional solution of
$(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi)=w$ in $\mathcal{O}$ with
$w\in\mathcal{L}^{\infty}_{loc}(\mathcal{O},\mathbb{R}^{m-1})$

  2. (2)

$\phi$ is locally $1/2$-Hölder continuous along the vertical components with
Hölder’s constant $C_{h}>0$

Then $\hat{\phi}$ is locally intrinsic Lipschitz function.

The proof of Theorem 7.2 relies on a preliminary result about the
distributional solutions of the problem
$(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi)=w$ in $\mathcal{O}$.

It is convenient to introduce the following notation: for $j=2,\dots,m$ and
$\hat{x}_{j}=(\hat{x}_{2},\dots,\hat{x}_{j-1},\hat{x}_{j+1},\dots,\hat{x}_{m})\in\mathbb{R}^{m-2}$
fixed, we denote by
$(t,\hat{x}_{j}):=(\hat{x}_{2},\dots,\hat{x}_{j-1},t,\hat{x}_{j+1},\dots,\hat{x}_{m})$

######  Lemma 7.3.

Under the same assumptions of Theorem 7.2, for all $j=2,\dots,m$ there is
$C=C(\|w_{j}\|_{\mathcal{L}^{\infty}(\mathcal{O},\mathbb{R})},C_{h},b_{j1}^{(1)},\dots,b_{j1}^{(n)})>0$
such that $\phi$ is a C-Lipschitz map along any characteristic line
$\gamma_{j}=(\gamma_{j1},\dots,\gamma_{jn}):[-T,T]\to\mathbb{R}^{n}$
satisfying

$\dot{\gamma}_{js}(t)=b_{j1}^{(s)}\phi(t,\hat{x}_{j},\gamma_{j}(t))+\frac{1}{2}\sum_{l=2}^{m}b^{(s)}_{jl}(\hat{x}_{j})_{l}\quad\mbox{
for all }s=1,\dots,n$

with $t\in[-T,T]$ and $\hat{x}_{j}\in\mathbb{R}^{m-2}$ fixed.

###### Proof.

Fix $j=2,\dots,m$ and $\hat{x}_{j}\in\mathbb{R}^{m-2}$.

Let $t_{1},t_{2}\in(-T,T)$ with $t_{1}<t_{2}$ and let $\epsilon>0$. In the
similar way in [9] and in [5], Lemma 5.1 we consider the functions
$\xi_{1}(x,y_{1}),\dots,\xi_{n}(x,y_{n})$ and $h(x)$ defined for small
$\delta>0$ by

$\xi_{s}(x,y_{s})=\left\\{\begin{array}[]{l}\quad\quad 0\hskip
91.04872pt-T<x<T,\,\,-r_{s}<y_{s}\leq\gamma_{s}(x)-\epsilon-\delta\\\
\frac{1}{\delta}(y_{s}-\gamma_{s}(x)+\epsilon+\delta)\hskip
18.49411pt-T<x<T,\,\,\gamma_{s}(x)-\epsilon-\delta<y_{s}\leq\gamma_{s}(x)-\epsilon\\\
\quad\quad 1\hskip
91.04872pt-T<x<T,\,\,\gamma_{s}(x)-\epsilon<y_{s}\leq\gamma_{s}(x)\\\
\frac{1}{\delta}(-y_{s}+\gamma_{s}(x)+\delta)\hskip
28.45274pt-T<x<T,\,\,\gamma_{s}(x)<y_{s}\leq\gamma_{s}(x)+\delta\\\ \quad\quad
0\hskip
91.04872pt-T<x<T,\,\,\gamma_{s}(x)+\delta<y_{s}<r_{s}\end{array}\right.$

for $s=1,\dots,n$ and

$h(x)=\left\\{\begin{array}[]{l}\quad\quad 0\hskip 71.13188pt-T<x\leq
t_{1}-\delta\\\ \frac{1}{\delta}(x-t_{1}+\delta)\hskip
36.98866ptt_{1}-\delta<x\leq t_{1}\\\ \quad\quad 1\hskip 72.55461ptt_{1}<x\leq
t_{2}\\\ \frac{1}{\delta}(-x+t_{2}+\delta)\hskip 29.87547ptt_{2}<x\leq
t_{2}+\delta\\\ \quad\quad 0\hskip
73.97733ptt_{2}+\delta<x<T\end{array}\right.$

If $\phi$ is distributional solution of
$(D_{2}^{\phi}\phi,\dots,D_{m}^{\phi}\phi)=w$ in $\mathcal{O}$, for all test
function $\varphi$ we have

$\displaystyle\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\dots\int_{-r_{n}}^{r_{n}}(\phi(x,\hat{x}_{j},y)\varphi_{x}(x,y)$
$\displaystyle+f_{1}(\phi(x,\hat{x}_{j},y))\varphi_{y_{1}}(x,y)+\dots+f_{n}(\phi(x,y))\varphi_{y_{n}}(x,y)$
$\displaystyle+w_{j}(x,\hat{x}_{j},y)\varphi(x,y))\,dxdy=0$

where

(59) $\displaystyle f_{s}(\phi(x_{j},\hat{x}_{j},y))$
$\displaystyle=\frac{1}{2}\biggl{(}b_{j1}^{(s)}\phi^{2}(x_{j},\hat{x}_{j},y)+\phi(x_{j},\hat{x}_{j},y)\sum_{l=2}^{m}b_{jl}^{(s)}(\hat{x}_{j})_{l}\biggl{)}$

for all $(x_{j},\hat{x}_{j},y)\in\mathcal{O}$ and $s=1,\dots,n$.

We choose the test function
$\varphi(x,y_{1},\dots,y_{n})=\xi_{1}(x,y_{1})\dots\xi_{n}(x,y_{n})h(x)$ and
compute the limit $\delta\to 0^{+}$:

$\displaystyle\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\dots\int_{-r_{n}}^{r_{n}}(\phi\varphi_{x}+f_{1}(\phi)(\xi_{1})_{y_{1}}\xi_{2}\dots\xi_{n}h+\dots+f_{n}(\phi)\xi_{1}\dots(\xi_{n})_{y_{n}}h+w_{j}\xi_{1}\dots\xi_{n}h)\,dxdy$

For simplicity we consider the case $n=2$. We compute the following limit:

$\displaystyle\lim_{\delta\to 0^{+}}$
$\displaystyle\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}\biggl{(}\phi\biggl{(}(\xi_{1})_{x}\xi_{2}h+\xi_{1}(\xi_{2})_{x}h+\xi_{1}\xi_{2}h_{x}\biggl{)}+f_{1}(\phi)(\xi_{1})_{y_{1}}\xi_{2}h$
$\displaystyle\quad+f_{2}(\phi)\xi_{1}(\xi_{2})_{y_{2}}h+w_{j}\xi_{1}\xi_{2}h\biggl{)}\,dxdy_{1}dy_{2}=:\lim_{\delta\to
0^{+}}I_{1}+I_{2}+I_{3}+I_{4}+I_{5}+I_{6}.$

Then

$\displaystyle\lim_{\delta\to 0^{+}}I_{1}$ $\displaystyle=\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}\phi(\xi_{1})_{x}\xi_{2}h\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{a}^{b}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}\dot{\gamma}_{1}(x)\Bigl{(}\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})\Bigl{)}\,dxdy_{2}$
$\displaystyle\lim_{\delta\to 0^{+}}I_{2}$ $\displaystyle=\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}\phi\xi_{1}(\xi_{2})_{x}h\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\dot{\gamma}_{2}(x)\Bigl{(}\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x))-\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)\Bigl{)}\,dxdy_{1}$
$\displaystyle\lim_{\delta\to 0^{+}}I_{3}$ $\displaystyle=\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}\phi\xi_{1}\xi_{2}h_{x}\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{\gamma_{1}(t_{1})-\epsilon}^{\gamma_{1}(t_{1})}\int_{\gamma_{2}(t_{1})-\epsilon}^{\gamma_{2}(t_{1})}\phi(t_{1},\hat{x}_{j},y_{1},y_{2})\,dy_{1}dy_{2}-\int_{\gamma_{1}(t_{2})-\epsilon}^{\gamma_{1}(t_{2})}\int_{\gamma_{2}(t_{2})-\epsilon}^{\gamma_{2}(t_{2})}\phi(t_{2},\hat{x}_{j},y_{1},y_{2})\,dy_{1}dy_{2}$
$\displaystyle\lim_{\delta\to 0^{+}}I_{4}$ $\displaystyle=\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}f_{1}(\phi)(\xi_{1})_{y_{1}}\xi_{2}h\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{t_{1}}^{t_{2}}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}\Bigl{(}f_{1}(\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2}))-f_{1}(\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2}))\Bigl{)}\,dxdy_{2}$
$\displaystyle\lim_{\delta\to 0^{+}}I_{5}$ $\displaystyle=\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}f_{2}(\phi)\xi_{1}(\xi_{2})_{y_{2}}h\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\Bigl{(}f_{2}(\phi(x,y_{1},\gamma_{2}(x)-\epsilon))-f_{2}(\phi(x,y_{1},\gamma_{2}(x)))\Bigl{)}\,dxdy_{1}$
$\displaystyle\lim_{\delta\to 0^{+}}I_{6}$ $\displaystyle=\lim_{\delta\to
0^{+}}\int_{-T}^{T}\int_{-r_{1}}^{r_{1}}\int_{-r_{2}}^{r_{2}}w_{j}\xi_{1}\xi_{2}h\,dxdy_{1}dy_{2}=\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}w_{j}\,dxdy_{1}dy_{2}$

Hence we obtain

(60)
$\displaystyle\int_{\gamma_{1}(t_{2})-\epsilon}^{\gamma_{1}(t_{2})}\int_{\gamma_{2}(t_{2})-\epsilon}^{\gamma_{2}(t_{2})}\phi(t_{2},\hat{x}_{j},y_{1},y_{2})\,dy_{1}dy_{2}-\int_{\gamma_{1}(t_{1})-\epsilon}^{\gamma_{1}(t_{1})}\int_{\gamma_{2}(t_{1})-\epsilon}^{\gamma_{2}(t_{1})}\phi(t_{1},\hat{x}_{j},y_{1},y_{2})\,dy_{1}dy_{2}$
$\displaystyle\quad-\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}w_{j}(x,\hat{x}_{j},y_{1},y_{2})\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{t_{1}}^{t_{2}}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}f_{1}(\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2}))-f_{1}(\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2}))$
$\displaystyle\quad+\dot{\gamma}_{1}(x)\Bigl{(}\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})\Bigl{)}\,dxdy_{2}$
$\displaystyle\quad+\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}f_{2}(\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon))-f_{2}(\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)))$
$\displaystyle\quad+\dot{\gamma}_{2}(x)\Bigl{(}\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x))-\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)\Bigl{)}\,dxdy_{1}$

From (59) with $n=2$ we deduce that

(61)
$\displaystyle\int_{\gamma_{1}(t_{2})-\epsilon}^{\gamma_{1}(t_{2})}\int_{\gamma_{2}(t_{2})-\epsilon}^{\gamma_{2}(t_{2})}\phi(t_{2},\hat{x}_{j},y_{1},y_{2})\,dy_{1}dy_{2}-\int_{\gamma_{1}(t_{1})-\epsilon}^{\gamma_{1}(t_{1})}\int_{\gamma_{2}(t_{1})-\epsilon}^{\gamma_{2}(t_{1})}\phi(t_{1},\hat{x}_{j},y_{1},y_{2})\,dy_{1}dy_{2}$
$\displaystyle\quad-\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}w_{j}(x,\hat{x}_{j},y_{1},y_{2})\,dxdy_{1}dy_{2}$
$\displaystyle=\int_{t_{1}}^{t_{2}}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}\frac{1}{2}b_{j1}^{(1)}\biggl{(}\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})\biggl{)}\biggl{(}\phi(x,\hat{x}_{j},\hat{x}_{j},\gamma_{1}(x),y_{2})+$
$\displaystyle+\phi(x,\hat{x}_{j},\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})-2\phi(x,\hat{x}_{j},\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))\biggl{)}\,dxdy_{2}+\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\frac{1}{2}b_{j1}^{(2)}\biggl{(}\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)$
$\displaystyle-\phi(x,y_{1},\gamma_{2}(x))\biggl{)}\biggl{(}\phi(x,y_{1},\gamma_{2}(x))+\phi(x,y_{1},\gamma_{2}(x)-\epsilon)-2\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))\biggl{)}\,dxdy_{1}$

If $b_{j1}^{(1)},b_{j1}^{(2)}=0$ the integral (61) is zero. On the other hand
if $b_{j1}^{(1)},b_{j1}^{(2)}\neq 0$, because $\phi$ is locally $1/2$-Hölder
continuous along the vertical components with Hölder’s constant $C_{h}>0$

(62)
$\displaystyle\frac{|\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})|}{\sqrt{\epsilon}}\leq
C_{h}$
$\displaystyle\frac{|\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))|}{\sqrt{\epsilon}}\leq
C_{h}$
$\displaystyle\frac{|\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))|}{\sqrt{2\epsilon}}\leq
C_{h}$

for $x\in[t_{1},t_{2}]$ and $y_{2}\in(\gamma_{2}(x)-\epsilon,\gamma_{2}(x))$.
Hence dividing by $\epsilon^{2}$ in $\eqref{lipottobre}$ and getting to the
limit to $\epsilon\to 0$, we obtain

(63) $\displaystyle\lim_{\epsilon\to
0}\frac{1}{\epsilon^{2}}\int_{t_{1}}^{t_{2}}\int_{\gamma_{2}(x)-\epsilon}^{\gamma_{2}(x)}\frac{1}{2}b_{j1}^{(1)}\biggl{(}\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})-\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})\biggl{)}\biggl{(}\phi(x,\hat{x}_{j},\gamma_{1}(x),y_{2})+$
$\displaystyle\quad+\phi(x,\hat{x}_{j},\gamma_{1}(x)-\epsilon,y_{2})-2\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))\biggl{)}\,dxdy_{2}$
$\displaystyle\leq|b_{j1}^{(1)}|\frac{(1+\sqrt{2})C_{h}^{2}}{2}(t_{2}-t_{1}).$

Moreover

(64) $\displaystyle\lim_{\epsilon\to
0}\frac{1}{\epsilon^{2}}\int_{t_{1}}^{t_{2}}\int_{\gamma_{1}(x)-\epsilon}^{\gamma_{1}(x)}\frac{1}{2}b_{j1}^{(2)}\biggl{(}\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)-\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x))\biggl{)}\biggl{(}\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x))+$
$\displaystyle\quad+\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)-2\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))\biggl{)}\,dxdy_{1}$
$\displaystyle\leq|b_{j1}^{(2)}|\frac{(1+\sqrt{2})C_{h}^{2}}{2}(t_{2}-t_{1}),$

where in the last inequality we have used the fact

(65)
$\displaystyle\frac{|\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)-\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x))|}{\sqrt{\epsilon}}\leq
C_{h}$
$\displaystyle\frac{|\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x))-\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))|}{\sqrt{\epsilon}}\leq
C_{h}$
$\displaystyle\frac{|\phi(x,\hat{x}_{j},y_{1},\gamma_{2}(x)-\epsilon)-\phi(x,\hat{x}_{j},\gamma_{1}(x),\gamma_{2}(x))|}{\sqrt{2\epsilon}}\leq
C_{h}$

for $x\in[t_{1},t_{2}]$ and $y_{1}\in(\gamma_{1}(x)-\epsilon,\gamma_{1}(x))$.

Now putting together $\eqref{dise01}$, $\eqref{dise02}$ and $\eqref{dise03}$,
dividing by $\epsilon^{2}$ and getting to the limit to $\epsilon\to 0$, we
obtain

(66)
$\phi(t_{2},\hat{x}_{j},\gamma_{j}(t_{2}))-\phi(t_{1},\hat{x}_{j},\gamma_{j}(t_{1}))\leq\biggl{(}\|w_{j}\|_{\mathcal{L}^{\infty}(\mathcal{O},\mathbb{R})}+\frac{(1+\sqrt{2})C_{h}^{2}}{2}\Big{(}|b_{j1}^{(1)}|+|b_{j1}^{(2)}|\Big{)}\biggl{)}(t_{2}-t_{1})$

for $t_{1},t_{2}\in(-T,T)$ with $t_{1}<t_{2}$.

In the similar way, using again (62) and (65) we have

(67)
$\phi(t_{2},\hat{x}_{j},\gamma_{j}(t_{2}))-\phi(t_{1},\hat{x}_{j},\gamma_{j}(t_{1}))\geq-\biggl{(}\|w_{j}\|_{\mathcal{L}^{\infty}(\mathcal{O},\mathbb{R})}+\frac{(1+\sqrt{2})C_{h}^{2}}{2}\Big{(}|b_{j1}^{(1)}|+|b_{j1}^{(2)}|\Big{)}\biggl{)}(t_{2}-t_{1})$

for $t_{1},t_{2}\in(-T,T)$ with $t_{1}<t_{2}$.

Hence combining (66) and (67) we get

$|\phi(t_{2},\hat{x}_{j},\gamma_{j}(t_{2}))-\phi(t_{1},\hat{x}_{j},\gamma_{j}(t_{1}))|\leq\biggl{(}\|w_{j}\|_{\mathcal{L}^{\infty}(\mathcal{O},\mathbb{R})}+\frac{(1+\sqrt{2})C_{h}^{2}}{2}\Big{(}|b_{j1}^{(1)}|+|b_{j1}^{(2)}|\Big{)}\biggl{)}(t_{2}-t_{1}),$

for $t_{1},t_{2}\in(-T,T)$ with $t_{1}<t_{2}$, i.e. $\phi$ is a Lipschitz map
along any characteristic line $\gamma_{j}$ as desired.

∎

######  Example 7.1.

We consider free step 2 group $\mathbb{F}_{m,2}$. As we said in Remark 3.1,
the composition law $\eqref{1}$ is given by
$\mathcal{B}^{(s)}\equiv\mathcal{B}^{(l,h)}$ where $1\leq h<l\leq m$ and
$\mathcal{B}^{(l,h)}$ has entries $-1$ in position $(l,h)$, $1$ in position
$(h,l)$ and $0$ everywhere else.

For each $j=2,\dots,m$ and $\hat{x}_{j}\in\mathbb{R}^{m-2}$ fixed the
characteristic line
$\gamma_{j}=(\gamma_{j1},\dots,\gamma_{jn}):[-\delta,\delta]\to\mathbb{R}^{n}$
has the following form:

$\displaystyle\dot{\gamma}_{js_{1}}(t)$
$\displaystyle=b_{j1}^{(s_{1})}\phi(t,\hat{x}_{j},\gamma_{j}(t))=-\phi(t,\hat{x}_{j},\gamma_{j}(t))$
$\displaystyle\dot{\gamma}_{js_{2}}(t)$
$\displaystyle=\frac{1}{2}b^{(s_{2})}_{jh_{2}}(\hat{x}_{j})_{h_{2}}=-\frac{1}{2}(\hat{x}_{j})_{h_{2}}$
$\displaystyle\dot{\gamma}_{js_{3}}(t)$ $\displaystyle=0$

where $\mathcal{B}^{(s_{i})}\equiv\mathcal{B}^{(l_{i},h_{i})}$ is such that

  * •

$s_{1}$ is the unique index $(l_{1},h_{1})=(j,1)$.

  * •

$s_{2}$ is such that $h_{2}<l_{2}=j$ and $h_{2}\neq 1$.

  * •

$s_{3}$ represents the other cases.

Now we are able to show the proof of Theorem 7.2.

###### Proof.

Fix $b\in\mathcal{O}$. According to Remark 2.15, we would like to prove that
there exists $\hat{C}>0$ such that

(68)
$|\phi(a^{\prime})-\phi(a)|\leq\hat{C}\|\hat{\phi}(i(a))^{-1}i(a)^{-1}i(a^{\prime})\hat{\phi}(i(a))\|\qquad\mbox{for
all }\,a,a^{\prime}\in\mathcal{U}(b,\delta)$

with $\delta>0$. Let $a=(x,y),a^{\prime}=(x^{\prime},y^{\prime})$ points of
$\mathcal{O}$ be sufficiently close to $b$, and let $\bar{D}_{j}$ be the
vector fields given by $\bar{D}_{j}=(x^{\prime}_{j}-x_{j})D^{\phi}_{j}$ for
$j\in\\{2,\dots,m\\}$. We define

$\displaystyle a_{1}$ $\displaystyle:=a$ $\displaystyle a_{2}$
$\displaystyle:=\exp(\bar{D}_{2})(a_{1})$ $\displaystyle a_{3}$
$\displaystyle:=\exp(\bar{D}_{3})(a_{2})$ $\displaystyle\vdots$ $\displaystyle
a_{m}$ $\displaystyle:=\exp(\bar{D}_{m})(a_{m-1}).$

More precisely for $j\in\\{2,\dots,m\\}$

$a_{j}=(x^{\prime}_{2},\dots,x^{\prime}_{j},x_{j+1},\dots,x_{m},y^{a_{j}})$

with

$\displaystyle y^{a_{j}}_{s}$
$\displaystyle=y_{s}+\sum_{l=2}^{j}\left(b_{l1}^{(s)}\int_{0}^{x^{\prime}_{l}-x_{l}}\phi\Big{(}\exp(rD^{\phi}_{l}(a_{l-1}))\Big{)}\,dr+\frac{1}{2}(x^{\prime}_{l}-x_{l})\Big{(}\sum_{h=2}^{l}x^{\prime}_{h}b_{lh}^{(s)}+\sum_{h=l+1}^{m}x_{h}b_{lh}^{(s)}\Big{)}\right)$
$\displaystyle=y^{a_{j-1}}_{s}+\left(b_{j1}^{(s)}\int_{0}^{x^{\prime}_{j}-x_{j}}\phi\Big{(}\exp(rD^{\phi}_{j}(a_{j-1}))\Big{)}\,dr+\frac{1}{2}(x^{\prime}_{j}-x_{j})\Big{(}\sum_{h=2}^{j}x^{\prime}_{h}b_{jh}^{(s)}+\sum_{h=j+1}^{m}x_{h}b_{jh}^{(s)}\Big{)}\right)$

for $s=1,\dots,n$. By Lemma 7.3 $(-\delta,\delta)\ni
r\mapsto\phi\Big{(}\exp(rD^{\phi}_{j}(a_{j-1}))\Big{)}$ is Lipschitz for all
$j=2,\dots,m$ and so $a_{2},\dots,a_{m}$ are well defined if
$a,a^{\prime}\in\mathcal{U}(b,\delta)$ for a sufficiently small $\delta$.

We observe that

(69)
$|\phi(a^{\prime})-\phi(a)|\leq|\phi(a^{\prime})-\phi(a_{m})|+\sum_{l=2}^{m}|\phi(a_{l})-\phi(a_{l-1})|$

and by Lemma 7.3 there is $C>0$ such that

(70)
$\sum_{l=2}^{m}|\phi(a_{l})-\phi(a_{l-1})|\leq\sum_{l=2}^{m}C|x^{\prime}_{l}-x_{l}|_{\mathbb{R}^{m-1}}\leq
C\|\hat{\phi}(i(a))^{-1}i(a)^{-1}i(a^{\prime})\hat{\phi}(i(a))\|.$

Hence in order to establish $\eqref{lipsch}$ we show that there exists
$C_{1}>0$ such that

(71) $|\phi(a^{\prime})-\phi(a_{m})|\leq
C_{1}\|\hat{\phi}(i(a))^{-1}i(a)^{-1}i(a^{\prime})\hat{\phi}(i(a))\|.$

Recalling that $a_{m}=(x^{\prime},y^{a_{m}})$ and using $\phi$ is locally
$1/2$-Hölder continuous along the vertical components with Hölder’s constant
$C_{h}>0$, it follows

$|\phi(a^{\prime})-\phi(a_{m})|\leq
C_{h}|y^{\prime}-y^{a_{m}}|^{1/2}_{\mathbb{R}^{n}}$

Moreover arguing as in the proof of the implication $4.\Rightarrow 2.$ in
Theorem 5.7 in [10], we have

$\displaystyle|y^{\prime}-y^{a_{m}}|_{\mathbb{R}^{n}}$
$\displaystyle\leq\sum_{s=1}^{n}\biggl{|}y^{\prime}_{s}-y_{s}+\sum_{l=2}^{m}\biggl{(}b_{1l}^{(s)}\int_{0}^{x^{\prime}_{l}-x_{l}}\phi\Big{(}\exp(rD^{\phi}_{l}(a_{l-1}))\Big{)}\,dr+$
$\displaystyle\quad-\frac{1}{2}(x^{\prime}_{l}-x_{l})\Big{(}\sum_{i=2}^{l}x^{\prime}_{i}b_{li}^{(s)}+\sum_{i=l+1}^{m}x_{i}b_{li}^{(s)}\Big{)}\biggl{)}\biggl{|}$
$\displaystyle\leq\sum_{s=1}^{n}\biggl{|}y^{\prime}_{s}-y_{s}+\phi(a)\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})b_{1l}^{(s)}-\frac{1}{2}\langle\mathcal{B}^{(s)}x,x^{\prime}-x\rangle\biggl{|}$
$\displaystyle\quad+\sum_{s=1}^{n}\biggl{|}-\frac{1}{2}\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})\Big{(}\sum_{i=2}^{l}x^{\prime}_{i}b_{li}^{(s)}+\sum_{i=l+1}^{m}x_{i}b_{li}^{(s)}\Big{)}+\frac{1}{2}\langle\mathcal{B}^{(s)}x,x^{\prime}-x\rangle\biggl{|}$
$\displaystyle\quad+\sum_{s=1}^{n}\biggl{|}-\phi(a)\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})b_{1l}^{(s)}+\sum_{l=2}^{m}b_{1l}^{(s)}\int_{0}^{x^{\prime}_{l}-x_{l}}\phi\Big{(}\exp(rD^{\phi}_{l}(a_{l-1}))\Big{)}\,dr\biggl{|}$
$\displaystyle\leq
c_{1}\|\hat{\phi}(i(a))^{-1}i(a)^{-1}i(a^{\prime})\hat{\phi}(i(a))\|^{2}+\frac{1}{2}n\mathcal{B}_{M}|x^{\prime}-x|^{2}_{\mathbb{R}^{m-1}}+$
$\displaystyle\quad+\sum_{s=1}^{n}\biggl{|}-\phi(a)\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})b_{1l}^{(s)}+\sum_{l=2}^{m}b_{1l}^{(s)}\int_{0}^{x^{\prime}_{l}-x_{l}}\phi\Big{(}\exp(rD^{\phi}_{l}(a_{l-1}))\Big{)}\,dr\biggl{|}$

where $c_{1}$ is given by (6) and
$\mathcal{B}_{M}=\max\\{b_{ij}^{(s)}\,|\,i,j=1,\dots,m\,,s=1,\dots,n\\}$. Note
that we have used

$\displaystyle\frac{1}{2}\langle\mathcal{B}^{(s)}x,x^{\prime}$
$\displaystyle-x\rangle-\frac{1}{2}\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})\Big{(}\sum_{i=2}^{l}x^{\prime}_{i}b_{li}^{(s)}+\sum_{i=l+1}^{m}x_{i}b_{li}^{(s)}\Big{)}$
$\displaystyle=-\frac{1}{2}\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})\Big{(}\sum_{i=2}^{l}x^{\prime}_{i}b_{li}^{(s)}+\sum_{i=l+1}^{m}x_{i}b_{li}^{(s)}-\sum_{i=2}^{m}x_{i}b_{li}^{(s)}\Big{)}$
$\displaystyle\leq\frac{1}{2}n\mathcal{B}_{M}|x^{\prime}-x|^{2}_{\mathbb{R}^{m-1}}.$

Finally, the last term

$\displaystyle\sum_{s=1}^{n}\biggl{|}-\phi(a)\sum_{l=2}^{m}(x^{\prime}_{l}-x_{l})b_{1l}^{(s)}+\sum_{l=2}^{m}b_{1l}^{(s)}\int_{0}^{x^{\prime}_{l}-x_{l}}\phi\Big{(}\exp(rD^{\phi}_{l}(B_{l-1}))\Big{)}\,dr\biggl{|}$
$\displaystyle\qquad\leq R_{1}(a,a^{\prime})+R_{2}(a,a^{\prime})$

where

$\displaystyle R_{1}(a,a^{\prime})$
$\displaystyle:=\sum_{s=1}^{n}\sum_{l=2}^{m}\biggl{|}b_{1l}^{(s)}\int_{0}^{x^{\prime}_{l}-x_{l}}\phi\Big{(}\exp(rD^{\phi}_{l}(a_{l-1}))\Big{)}\,dr-
b_{1l}^{(s)}\phi(a_{l-1})(x^{\prime}_{l}-x_{l})\biggl{|}$ $\displaystyle
R_{2}(a,a^{\prime})$
$\displaystyle:=\sum_{s=1}^{n}\biggl{|}\sum_{l=2}^{m}b_{1l}^{(s)}(x^{\prime}_{l}-x_{l})\Big{(}\phi(a_{l-1})-\phi(a)\Big{)}\biggl{|}$

We would like to show that there exist $C_{1},C_{2}>0$ such that

(72) $R_{1}(a,a^{\prime})\leq C_{1}|x^{\prime}-x|^{2}_{\mathbb{R}^{m-1}}$ (73)
$R_{2}(a,a^{\prime})\leq C_{2}|x^{\prime}-x|^{2}_{\mathbb{R}^{m-1}}$

for all $a,a^{\prime}\in\mathcal{U}(b,\delta)$ and consequently

$|y^{\prime}-y^{a_{m}}|_{\mathbb{R}^{n}}\leq
c_{1}\|\hat{\phi}(i(a))^{-1}i(a)^{-1}i(a^{\prime})\hat{\phi}(i(a))\|^{2}+\left(\frac{1}{2}n\mathcal{B}_{M}+C_{1}+C_{2}\right)|x^{\prime}-x|^{2}_{\mathbb{R}^{m-1}}$

Hence there is $C_{3}>0$ such that

(74) $|y^{\prime}-y^{a_{m}}|^{1/2}_{\mathbb{R}^{n}}\leq
C_{3}\|\hat{\phi}(i(a))^{-1}i(a)^{-1}i(a^{\prime})\hat{\phi}(i(a))\|,$

i.e. (71) is true. We start to consider $R_{1}(a,a^{\prime})$. Fix
$l=2,\dots,m$. For $t\in[-\delta,\delta]$ we define

$g_{l}(t):=\sum_{s=1}^{n}b_{1l}^{(s)}\biggl{(}\int_{0}^{t}\phi(\exp(rD^{\phi}_{l})(a_{l-1}))\,dr-t\phi(a_{l-1})\biggl{)}$

Observe that

$\displaystyle
b_{1l}^{(s)}\int_{0}^{t}\left(\phi(\exp(rD^{\phi}_{l})(a_{l-1}))-\phi(a_{l-1})\right)\,dr=O(t^{2})$

and so there is $C_{l}>0$ such that

$|g_{l}(t)|\leq C_{l}t^{2},\hskip 14.22636pt\forall t\in[-\delta,\delta]$

Hence set $t=x^{\prime}_{l}-x_{l}$ we get

$|g_{l}(x^{\prime}_{l}-x_{l})|\leq C_{l}(x^{\prime}_{l}-x_{l})^{2}$

and consequently $\eqref{r1}$ follows from

$\sum_{l=2}^{m}|g_{l}(x^{\prime}_{l}-x_{l})|\leq\sum_{l=2}^{m}C_{l}(x^{\prime}_{l}-x_{l})^{2}\leq
C_{1}|x^{\prime}-x|_{\mathbb{R}^{m-1}}^{2}.$

Now we consider $R_{2}(a,a^{\prime})$. Observe that

$\displaystyle\sum_{s=1}^{n}\biggl{|}\sum_{l=2}^{m}b_{1l}^{(s)}(x^{\prime}_{l}$
$\displaystyle-x_{l})\Big{(}\phi(a_{l-1})-\phi(a)\Big{)}\biggl{|}\leq
n\mathcal{B}_{M}\sum_{l=2}^{m}|x^{\prime}_{l}-x_{l}|\left|\phi(a_{l-1})-\phi(a)\right|$
$\displaystyle\leq
n\mathcal{B}_{M}\sum_{l=2}^{m}|x^{\prime}_{l}-x_{l}|\Big{(}\sum_{h=2}^{l-1}|\phi(a_{h})-\phi(a_{h-1})|\Big{)}$
$\displaystyle\leq
n\mathcal{B}_{M}\sum_{l=2}^{m}|x^{\prime}_{l}-x_{l}|\Big{(}\sum_{h=2}^{l-1}\Big{|}\int_{0}^{1}(\bar{D}_{h}\phi)(\exp(r\bar{D}_{h}(a_{h-1})))\,dr\Big{|}\Big{)}$
$\displaystyle\leq
n\mathcal{B}_{M}\sum_{l=2}^{m}|x^{\prime}_{l}-x_{l}|\Big{(}\sum_{h=2}^{l-1}\Big{|}(x^{\prime}_{h}-x_{h})\left(D^{\phi}_{h}\phi(a_{h-1})+o(1)\right)\Big{|}\Big{)}$
$\displaystyle\leq n\mathcal{B}_{M}C|x^{\prime}-x|^{2}_{\mathbb{R}^{m-1}}$

Then $\eqref{r2}$ follows with $C_{2}:=n\mathcal{B}_{M}C$ and (74) is true.

Finally putting together $\eqref{rimanente5}$, $\eqref{rimanente4}$ and
$\eqref{rimanenteok1}$, $\eqref{lipsch}$ holds and the proof is complete. ∎

##  8\. A characterization of $\mathbb{G}$-regular hypersurfaces

In this section we will prove a new characterization of $\mathbb{G}$-regular
hypersurfaces (see Definition 2.2). Here we consider the non linear first
order system

(75) $\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w$

where $w:\mathcal{O}\to\mathbb{R}^{m-1}$ is a given _continuous function_ and
not measurable as opposed to the previous sections.

In [10], we give some equivalent conditions about $\mathbb{G}$-regular
hypersurface inside the Carnot groups of step 2 (see Theorem 5.7 in [10]). In
this section we give another one. More precisely, we show that if a continuous
map $\phi:\mathcal{O}\subset\mathbb{R}^{m+n-1}\to\mathbb{R}$ is locally little
$1/2$-Hölder continuous along the vertical components, then $\phi$ is a
distributional solution in an open set $\mathcal{O}$ of the non linear first
order system $\eqref{sistema}$ if and only if
$\hat{\phi}:\hat{\mathcal{O}}\to\mathbb{V}$ is uniformly intrinsic
differentiable in $\hat{\mathcal{O}}$ and consequently its graph is a
$\mathbb{G}$-regular hypersurface. The main equivalence result is contained in
Theorem 8.2.

In order to state the equivalence result we have to be precise about the
meaning of being a solution of (75). To this aim we recall a notion of
generalized solutions of systems of this kind. These generalized solutions,
denoted _broad* solutions_ were introduced and studied for the system (75)
inside Heisenberg groups in [3, 6]. For a more complete bibliography we refer
to the bibliography in [3]. Our strategy will be to prove that each continuous
distributional solution of the system (75) is a broad* solution and vice versa
(see Proposition 8.1).

######  Definition 8.1.

Let $\mathcal{O}\subset\mathbb{R}^{m+n-1}$ be open and
$w:=\left(w_{2},\dots,w_{m}\right):\mathcal{O}\to\mathbb{R}^{m-1}$ a
continuous function. With the notations of Definition 3.1 we say that
$\phi\in\mathbb{C}(\mathcal{O},\mathbb{R})$ is a _broad* solution_ in
$\mathcal{O}$ of the system

$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w$

if for every $a_{0}\in\mathcal{O}$ there are $0<\delta_{2}<\delta_{1}$ and
$m-1$ maps $\exp_{a_{0}}(\cdot D^{\phi}_{j})(\cdot)$

$\begin{split}\exp_{a_{0}}(\cdot
D^{\phi}_{j})(\cdot):[-\delta_{2},\delta_{2}]&\times\mathcal{I}(a_{0},\delta_{2})\to\mathcal{I}(a_{0},\delta_{1})\\\
(t,&a)\mapsto\exp_{a_{0}}(tD^{\phi}_{j})(a)\end{split}$

for $j=2,\dots,m$, where
$\mathcal{I}(a_{0},\delta):=\mathcal{U}(a_{0},\delta)\cap\mathbb{W}$ and
$\mathcal{I}(a_{0},\delta_{1})\subset\mathcal{O}$. Moreover these maps, called
_exponential maps_ of the vector fields $D^{\phi}_{2},\dots,D^{\phi}_{m}$, are
required to have the following properties

$t\mapsto\gamma^{j}_{a}(t):=\exp_{a_{0}}(tD^{\phi}_{j})(B)\in\mathbb{C}^{1}([-\delta_{2},\delta_{2}],\mathbb{R}^{m+n-1})$

for all $a\in\mathcal{I}(a_{0},\delta_{2})$ and

$\left\\{\begin{array}[]{l}\dot{\gamma}^{j}_{a}=D^{\phi}_{j}\circ\gamma^{j}_{a}\\\
\gamma^{j}_{a}(0)=a\end{array}\right.$
$\phi(\gamma^{j}_{a}(t))-\phi(\gamma^{j}_{a}(0))=\int_{0}^{t}w_{j}(\gamma^{j}_{a}(r))\,dr$

once more for all $a\in\mathcal{I}(a_{0},\delta_{2})$.

######  Proposition 8.1.

Let $\mathbb{G}:=(\mathbb{R}^{m+n},\cdot,\delta_{\lambda})$ be a Carnot group
of step $2$ and $\mathbb{V}$, $\mathbb{W}$ the complementary subgroups defined
in (8). Let $\hat{\phi}:\hat{\mathcal{O}}\to\mathbb{V}$ be a continuous map,
where $\hat{\mathcal{O}}$ is an open subset of $\mathbb{W}$ and
$\phi:\mathcal{O}\to\mathbb{R}$ is the map associated to $\hat{\phi}$ as in
$\eqref{phipsi}$. Then the following conditions are equivalent:

  1. (1)

$\phi$ is locally little $1/2$-Hölder continuous along the vertical
components, i.e. that is $\phi\in\mathbb{C}(\mathcal{O},\mathbb{R})$ and for
all $\mathcal{O}^{\prime}\Subset\mathcal{O}$ and
$(x,y),(x,y^{\prime})\in\mathcal{O}^{\prime}$

$\lim_{r\to
0^{+}}\sup_{0<|y^{\prime}-y|_{\mathbb{R}^{n}}<r}\,\frac{|\phi(x,y^{\prime})-\phi(x,y)|}{|y^{\prime}-y|_{\mathbb{R}^{n}}^{1/2}}=0.$

and there exists $w\in\mathbb{C}(\mathcal{O},\mathbb{R}^{m-1})$ such that
$\phi$ is a broad* solution of

$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w,\quad\mbox{ in
}\mathcal{O}$

  2. (2)

$\phi$ is locally little $1/2$-Hölder continuous along the vertical components
and there exists $w\in\mathbb{C}(\mathcal{O},\mathbb{R}^{m-1})$ such that
$\phi$ is distributional solution of

$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w,\quad\mbox{ in
}\mathcal{O}$

###### Proof.

$\mathbf{(1)\implies(2)}$. By Theorem $5.7$ in [10], we know that there is a
family of functions $\phi_{\epsilon}\in\mathbb{C}^{1}(\mathcal{O},\mathbb{R})$
such that for all $\mathcal{O}^{\prime}\Subset\mathcal{O}$,

(76) $\phi_{\epsilon}\to\phi\quad and\quad
D_{j}^{{\phi_{\epsilon}}}{\phi_{\epsilon}}\to D_{j}^{{\phi}}{\phi}$

for $j=2,\dots,m$ uniformly on $\mathcal{O}^{\prime}$ as $\epsilon\to 0^{+}$.
Then for each $j=2,\dots,m$, $\epsilon>0$ and
$\zeta\in\mathbb{C}^{1}_{c}(\mathcal{O},\mathbb{R})$

$\int_{\mathcal{O}}\phi_{\epsilon}\left(\,X_{j}\zeta+\phi_{\epsilon}\sum_{s=1}^{n}b_{j1}^{(s)}Y_{s}\zeta\right)\,d\mathcal{L}^{m+n-1}=-\int_{\mathcal{O}}D^{\phi_{\epsilon}}_{j}\phi_{\epsilon}\,\zeta\,d\mathcal{L}^{m+n-1}$

Using (76) and getting to the limit for $\epsilon\to 0^{+}$ we have that

$\int_{\mathcal{O}}\phi\left(\,X_{j}\zeta+\phi\sum_{s=1}^{n}b_{j1}^{(s)}Y_{s}\zeta\right)\,d\mathcal{L}^{m+n-1}=-\int_{\mathcal{O}}D^{\phi}_{j}\phi\,\zeta\,d\mathcal{L}^{m+n-1}$

Hence $\phi$ is distributional solution of
$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w$ in $\mathcal{O}$.

$\mathbf{(2)\implies(1)}$. Let $a_{0}\in\mathcal{O}$, $\delta_{1}>0$ and
$\mathcal{I}:=\mathcal{U}(a_{0},\delta_{1})\cap\mathbb{W}$ and
$\mathcal{I}(a_{0},\delta_{1})\subset\mathcal{O}$. Denote

$K:=\sup_{(x,y)\in\mathcal{I}}\sum_{l=2}^{m}|x_{l}|,\quad
M:=\|\phi\|_{\mathcal{L}^{\infty}(\mathcal{I},\mathbb{R})},\quad\delta_{2}<\frac{\delta_{1}}{2+\frac{1}{2}K\mathcal{B}_{M}+M\mathcal{B}_{M}},$

where
$\mathcal{B}_{M}=\max\\{b_{jl}^{(s)}\,|\,l,j=1,\dots,m\,,s=1,\dots,n\\}$.
Peano’s Theorem yields that for all
$a=(x_{j},\hat{x}_{j},y)\in\overline{\mathcal{I}(a_{0},\delta_{2})}$ there
exists $\mathbb{C}^{1}$ function
$(\gamma_{j1},\dots,\gamma_{jn}):[-\delta_{2},\delta_{2}]\to\mathbb{R}^{n}$
such that

$\gamma^{j}_{a}(t)=(x_{j}+t,\hat{x}_{j},\gamma_{j1}(t),\dots,\gamma_{jn}(t))\in\mathcal{I}\quad\mbox{for
}t\in[-\delta_{2},\delta_{2}],$

and $\gamma_{js}(t)$ is a solution of the Cauchy problem

$\left\\{\begin{array}[]{l}\dot{\gamma}_{js}(t)=\frac{1}{2}\sum_{l=2}^{m}(\hat{x}_{j})_{l}b_{jl}^{(s)}+b_{j1}^{(s)}\phi(\gamma_{a}^{j}(t)),\quad\mbox{
for }t\in[-\delta_{2},\delta_{2}]\\\ \\\ \gamma_{js}(0)=y_{s},\\\
\end{array}\right.$

for $s=1,\dots,n$. It is clear that if $b_{j1}^{(s)}=0$, then
$\gamma_{js}(t)=y_{s}+\frac{1}{2}t\sum_{l=2}^{m}(\hat{x}_{j})_{l}b_{jl}^{(s)}$
for $t\in[-\delta_{2},\delta_{2}]$. On the other hand, if $b_{j1}^{(s)}\neq
0$, then the map $\phi(\gamma_{a}^{j}(\cdot))$ satisfies the following ODE

$\frac{d}{dt}\left(b_{j1}^{(s)}\phi(\gamma_{a}^{j}(t))\right)=b_{j1}^{(s)}w_{j}(\gamma_{a}^{j}(t))\\\
$

with $t\in[-\delta,\delta]$ for some $\delta>0$. Indeed, we can repeat
verbatim the proof of Lemma 7.3 and so we obtain (61) with $b_{j1}^{(s)}\neq
0$. Moreover because $\phi$ is locally little $1/2$-Hölder continuous along
the vertical components, upon dividing (61) by $\epsilon^{2}$ and getting to
the limit to $\epsilon\to 0$, we have that

$\displaystyle\phi(\gamma_{a}^{j}(t_{2}))-\phi(\gamma_{a}^{j}(t_{1}))-\int_{t_{1}}^{t_{2}}w_{j}(\gamma_{a}^{j}(t))\,dt=0$

for $t_{1},t_{2}\in[-\delta,\delta]$ with $\delta>0$ and $t_{1}<t_{2}$.

In particular $\phi(\gamma_{a}^{j}(\cdot))$ and $\dot{\gamma}_{js}(\cdot)$ are
$\mathbb{C}^{1}([-\delta_{3},\delta_{3}],\mathbb{R})$ for $s=1,\dots,n$ where
$\delta_{3}:=\min\\{\delta,\delta_{2}\\}$. Therefore the curve
$\gamma^{j}_{a}:[-\delta_{3},\delta_{3}]\to\mathcal{I}$ satisfies the
conditions of the Definition 8.1 for each
$a\in\mathcal{I}(a_{0},\delta_{3}):=\mathcal{U}(a_{0},\delta_{3})\cap\mathbb{W}$
and $\mathcal{I}(a_{0},\delta_{3})\subset\mathcal{O}$.

Then, for each $j=2,\dots,m$, $\exp_{a_{0}}(\cdot
D^{\phi}_{j})(\cdot):[-\delta_{3},\delta_{3}]\times\mathcal{I}(a_{0},\delta_{3})\to\mathcal{I}$
defined as $\exp_{a_{0}}(tD^{\phi}_{j})(a):=\gamma^{j}_{a}(t)$ is a family of
exponential maps at $a_{0}$ which we were looking for. This completes the
proof of the implication $(4)\implies(5)$.

∎

By Proposition 8.1 and Theorem $5.7$ in [10], it immediately follows

######  Theorem 8.2.

Let $\mathbb{G}:=(\mathbb{R}^{m+n},\cdot,\delta_{\lambda})$ be a Carnot group
of step $2$ and $\mathbb{V}$, $\mathbb{W}$ the complementary subgroups defined
in (8). Let $\hat{\phi}:\hat{\mathcal{O}}\to\mathbb{V}$ be a continuous map,
where $\hat{\mathcal{O}}$ is an open subset of $\mathbb{W}$ and
$\phi:\mathcal{O}\to\mathbb{R}$ is the map associated to $\hat{\phi}$ as in
$\eqref{phipsi}$. Then the following conditions are equivalent:

  1. (1)

$\mathrm{graph}\,(\hat{\phi})$ is a $\mathbb{G}$-regular hypersurface and for
all $a\in\mathrm{graph}\,(\hat{\phi})$ there is $r=r(a)>0$ and
$f\in\mathbb{C}_{\mathbb{G}}^{1}(\mathcal{U}(a,r)),\mathbb{R})$ with
$X_{1}f>0$ such that

$\mathrm{graph}\,(\hat{\phi})\cap\mathcal{U}(a,r)=\\{p:f(p)=0\\}$

  2. (2)

$\hat{\phi}$ is u.i.d. in $\hat{\mathcal{O}}$.

  3. (3)

$D_{j}^{\phi}\phi$ interpreted in distributional sense is a continuos function
in $\mathcal{O}$ and for $0<\epsilon<1$ there is a family of functions
$\phi_{\epsilon}\in\mathbb{C}^{1}(\mathcal{O},\mathbb{R})$ such that for all
$\mathcal{O}^{\prime}\Subset\mathcal{O}$,

$\phi_{\epsilon}\to\phi\quad and\quad
D_{j}^{{\phi_{\epsilon}}}{\phi_{\epsilon}}\to D_{j}^{{\phi}}{\phi}$

for $j=2,\dots,m$, uniformly on $\mathcal{O}^{\prime}$ as $\epsilon\to 0^{+}$.

  4. (4)

$\phi$ is locally little $1/2$-Hölder continuous along the vertical components
and there exists $w\in\mathbb{C}(\mathcal{O},\mathbb{R}^{m-1})$ such that
$\phi$ is a broad* solution of

$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w,\quad\mbox{ in
}\mathcal{O}$

  5. (5)

$\phi$ is locally little $1/2$-Hölder continuous along the vertical components
and there exists $w\in\mathbb{C}(\mathcal{O},\mathbb{R}^{m-1})$ such that
$\phi$ is distributional solution of

$\left(D^{\phi}_{2}\phi,\dots,D^{\phi}_{m}\phi\right)=w,\quad\mbox{ in
}\mathcal{O}$

######  Remark 8.3.

We recall that in Theorem $5.7$ in [10], we consider a locally little
$1/2$-Hölder continuous map but in the proof we just use that $\phi$ is
locally little $1/2$-Hölder continuous along the vertical components.

## References

  * [1] L.Ambrosio, _Some fine properties of sets of finite perimeter in Ahlfors regular metric measure spaces_ , Adv. in Math., $\mathbf{159}$, (2001), 51-67.
  * [2] L.Ambrosio, B.Kirchheim, _Rectifiable sets in metric and Banach spaces_ , Math. Ann. 318 (2000), no. 3, 527-555.
  * [3] L.Ambrosio, F.Serra Cassano, D.Vittone, _Intrinsic regular hypersurfaces in Heisenberg groups_ , J. Geom. Anal. $\mathbf{16}$, (2006), 187-232.
  * [4] G.Arena, R.Serapioni, _Intrinsic regular submanifolds in Heisenberg groups are differentiable graphs_ , Calc. Var. Partial Differential Equations $\mathbf{35}$, (2009), no. 4, 517-536.
  * [5] F.Bigolin, L.Caravenna, F.Serra Cassano, _Intrinsic Lipschitz graphs in Heisenberg groups and continuous solutions of a balance equation_ , Ann. I.H. Poincaré -AN (2014).
  * [6] F.Bigolin, F.Serra Cassano, _Intrinsic regular graphs in Heisenberg groups vs. weak solutions of non linear first-order PDEs_ , Adv. Calc. Var. 3, (2010), 69-97.
  * [7] A.Bonfiglioli, E.Lanconelli, F.Uguzzoni, _Stratified Lie Groups and Potential Theory for their Sub-Laplacians_ , Springer Monographs in Mathematics, Springer-Verlag Berlin Heidelberg, New York (2007).
  * [8] G.Citti, M.Manfredini, A.Pinamonti, F.Serra Cassano, _Smooth approximation for the intrinsic Lipschitz functions in the Heisenberg group_ , Calc. Var. Partial Differ. Equ. 49 (2014), 1279-1308.
  * [9] C.M. Dafermos, Continuous solutions for balance laws, Ric. Mat. 55 (2006), 79-91.
  * [10] D.Di Donato, _Intrinsic differentiability and Intrinsic Regular Surfaces in Carnot groups_ , arXiv:1811.05457, preprint.
  * [11] D.Di Donato, _Intrinsic differentiability and Intrinsic Regular Surfaces in Carnot groups_ , PhD Thesis, University of Trento, (2017).
  * [12] L.C.Evans, R.F.Gariepy, _Measure Theory and Fine Proprieties of functions_ , CRC Press, BocaRaton, (1992).
  * [13] H.Federer, _Geometric Measure Theory_ , Springer, (1969).
  * [14] G.B. Folland, E.M. Stein, _Hardy Spaces on Homogeneous Groups_ , Princeton University Press, (1982).
  * [15] B.Franchi, M.Marchi, R.Serapioni, _Differentiability and Approximate Differentiability for intrinsic Lipschitz functions in Carnot groups and a Rademarcher Theorem_ , Anal. Geom. Metr. Spaces $\mathbf{2}$, (2014), 258-281.
  * [16] B.Franchi, R.Serapioni, _Intrinsic Lipschitz Graph within Carnot groups_ , The Journal of Geometric Analysis, Vol. 26 (3) (2016), 1946-1994.
  * [17] B.Franchi, R.Serapioni, F.Serra Cassano, _Area Formula for Centered Hausdorff Measures in Metric Spaces_ , Nonlinear Analysis, (2015), 218-233.
  * [18] B.Franchi, R.Serapioni, F.Serra Cassano, _Intrinsic Lipschitz graphs in Heisenberg groups_ , J. Nonlinear Convex Anal. 7 (2006), no. 3, 423-441.
  * [19] B.Franchi, R.Serapioni, F.Serra Cassano, _Rectificability and perimeter in the Heisenberg group_ , Math. Ann. $\mathbf{321}$ (2001), 479-531.
  * [20] B.Franchi, R.Serapioni, F.Serra Cassano, _Regular hypersurfaces, Intrinsic Perimeter and Implicit Function Theorem in Carnot Groups_ , Comm. Anal. Geom. $\mathbf{11}$ (2003), no. 5, 909-944.
  * [21] B.Franchi, R.Serapioni, F.Serra Cassano, _Regular Submanifolds, Graphs and area formula in Heisenberg Groups_ , Advances in Math. $\mathbf{211}$ (2007), no. 1, 152-203.
  * [22] B.Franchi, R.Serapioni, F.Serra Cassano, _On the Structure of Finite Perimeter Sets in Step 2 Carnot Groups_ , The Journal of Geometric Analysis $\mathbf{13}$ (2003), no. 3, 421-466.
  * [23] N.Garofalo, D.M.Nhieu, _Isoperimetric and Sobolev inequalities for Carnot-Carathéodory spaces and the existence of minimal surfaces_ , Comm. Pure Appl. Math. $\mathbf{49}$ (1996), 1081-1144.
  * [24] M.Gromov, _Carnot-Carathéodory spaces seen from within, in Sub Riemannian Geometry_ , Progress in Mathematics, 144. ed. by A.Bellaïche and J.Risler, Birkhauser Verlag, Basel (1996).
  * [25] B.Kirchheim, F.Serra Cassano, _Rectifiability and parametrization of intrinsic regular surfaces in the Heisenberg group_ , Ann. Scuola Norm. Sup. Pisa Cl. Sci. (5) $\mathbf{III}$, (2004), 871-896.
  * [26] E.Le Donne, _A Primer on Carnot Groups: Homogenous Groups, Carnot-Carathéodory Spaces, and Regularity of Their Isometries_ , Anal. Geom. Metr. Spaces, (2017), 5: 116-137.
  * [27] V.Magnani, _Towards Differential Calculus in Stratified Groups_ , J. Aust. Math. Soc. $\mathbf{95}$, n.1, (2013), 76-128.
  * [28] P.Mattila, _Geometry of Sets and Measures in Euclidean Spaces_ , Cambridge University Press, Cambridge, (1995).
  * [29] J.Mitchell, _On Carnot-Carathéodory metrics_ , J. Differ. Geom. $\mathbf{21}$, (1985), 35-45.
  * [30] R.Monti, D.Vittone, _Sets with finite H-perimeter and controlled normal_ , Math. Z. 270 (2012), no. 1-2, 351-367.
  * [31] P.Pansu, _Métriques de Carnot-Carathéodory et quasiisométries des espaces symétriques de rang un_ , Ann. of Math., $\mathbf{129}$, (1989), 1-60.
  * [32] R.Serapioni, _Intrinsic difference quotients_ , Harmonic Analysis, Partial Differential Equations and Applications, Birkhäuser (2017), 165-192.
  * [33] F.Serra Cassano, _Geometry, Analysis and Dynamics on sub-Riemannian Manifolds_ , Vol. I, European Mathematical Society (2016).
  * [34] F.Serra Cassano, D.Vittone, _Graphs of bounded variation, existence and local boundedness of non-parametric minimal surfaces in the Heisenberg group_ , Adv. Calc. Var. 7 (2014), no. 4, 409-492.
  * [35] A.Visintin, _Strong Convergence results related to strict convexity_ , Comm. in PDE, 9 (1984), 439-466.
  * [36] D.Vittone, _Submanifolds in Carnot Groups_ , PhD Thesis, Scuola Normale Superiore, Edizioni della Normale, Pisa, (2008).