# ECRECer: Enzyme Commission Number Recommendation and Benchmarking based on
Multiagent Dual-core Learning ††thanks: Citation: Authors. Title. Pages….
DOI:000000/11111.

Zhenkun Shi, Qianqian Yuan, Ruoyu Wang, Hoaran Li, Xiaoping Liao🖂, Hongwu Ma🖂
Biodesign Center, Key Laboratory of Systems Microbial Technology
Tianjin Institute of Industrial Biotechnology, Chinese Academy of Sciences,
300308, Tianjin, China
National Technology Innovation Center of Synthetic Biology, 300308, Tianjin,
China
{Zhenkun.Shi, yuan_qq, wangry, lihr, liao_xp<EMAIL_ADDRESS>

###### Abstract

Enzyme Commission (EC) numbers, which associate a protein sequence with the
biochemical reactions it catalyzes, are essential for the accurate
understanding of enzyme functions and cellular metabolism. Many ab-initio
computational approaches were proposed to predict EC numbers for given input
sequences directly. However, the prediction performance (accuracy, recall,
precision), usability, and efficiency of existing methods still have much room
to be improved. Here, we report ECRECer, a cloud platform for accurately
predicting EC numbers based on novel deep learning techniques. To build
ECRECer, we evaluate different protein representation methods and adopt a
protein language model for protein sequence embedding. After embedding, we
propose a multi-agent hierarchy deep learning-based framework to learn the
proposed tasks in a multi-task manner. Specifically, we used an extreme multi-
label classifier to perform the EC prediction and employed a greedy strategy
to integrate and fine-tune the final model. Comparative analyses against four
representative methods demonstrate that ECRECer delivers the highest
performance, which improves accuracy and F1 score by 70% and 20% over the
state-of-the-the-art, respectively. With ECRECer, we can annotate numerous
enzymes in the Swiss-Prot database with incomplete EC numbers to their full
fourth level. Take UniPort protein "A0A0U5GJ41" as an example (1.14.-.-),
ECRECer annotated it with "1.14.11.38", which supported by further protein
structure analysis based on AlphaFold2. Finally, we established a webserver
(https://ecrecer.biodesign.ac.cn) and provided an offline bundle to improve
usability.

_K_ eywords EC prediction, protein language model, extreme multi-label
classification, deep learning

##  1 Introduction

With the widespread adoption of high-throughput methods and high-quality
infrastructure in biotechnology and bioindustry, the speed of new protein
discovery has increased dramatically. However, this was not followed by a
concomitant increase in the speed of protein annotation (see Supplemental, SI
Appendix FIGURES, Fig. S2). For example, 5 241 146 sequences were added to
TrEMBL in the UniProt database [1] in the single month of March 2021, while
only 521 sequences were reviewed and added to Swiss-Prot in the same period
(see Supplemental, SI Appendix FIGURES, Fig. S2). Such a slow speed of protein
annotation considerably restricts related research and industrial
applications.

Among the multiple and complex protein annotation tasks, one of the crucial
steps is enzyme function annotation [2, 3]. Annotations of enzyme function
provide critical starting points for generating and testing biological
hypotheses [3]. Current functional annotations of enzymes describe the
biochemistry or process by assigning an Enzyme Commission (EC) number. This is
a four-part code associated with a recommended name for the corresponding
enzyme-catalyzed reaction that describes the enzyme class, the chemical bond
acted on, the reaction, and the substrates [4]. Thus, the primary task of
enzyme annotation is to assign an EC number to a given protein sequence.
However, as the uncertainty of the assignments for uncharacterized protein
sequences is high and biochemical data are relatively sparse, both the speed
and the quality of enzyme annotation are considerably restricted.

To achieve improved, rapid and intelligent functional annotation,
computational methods were introduced to assign or predict EC numbers. The
simplest and most commonly used method is multiple sequence alignment (MSA)
[5], which can yield an appropriate annotation by using similar sequences.
Based on this approach, researchers have developed most major EC databases and
profile-based methods for the functional annotation of enzymes [6, 7, 8, 9].
However, these methods cannot perform annotations for novel proteins with no
similar sequences, which is generally the case for newly discovered enzymes.
To overcome this restriction, researchers introduced machine learning methods,
such as SVM [10], KNN [11], and hidden Markov model [12] for the functional
annotation of enzymes. Although these methods can predict EC numbers even if
the given protein sequences have no similar references, the prediction speed
and precision are not ideal. Since deep learning has delivered powerful
results in many areas [13, 14, 15, 16], more researchers are trying to use
deep learning methods to predict EC numbers and significantly improve the
precision of functional annotation. However, deep learning methods are prone
to overfitting due to an unbalanced distribution of training datasets. In EC
number prediction, this leads to prediction results with high precision,
medium recall, and low accuracy.

Overall, there has been a steady improvement in computational methods for
enzyme annotation [9, 7, 17, 2], but several obstacles still exist that have
slowed the progress of computational enzyme function annotation. One of the
direct challenges is a lack of publicly available benchmark datasets to
evaluate the existing and newly proposed models, which makes it troublesome
for the end-user to choose the best method in their production scenario.
Another notable challenge is the lack of an efficient and universal protein
sequence embedding method. Thus, researchers have to spend large amounts of
time on handcrafted feature engineering to encode the sequence, such as
functional domain encoding [18] and position-specific scoring matrix encoding
[19], as encoding quality dramatically impacts the performance of downstream
applications [20]. The third challenge is the lack of an explicitly designed
method to deal with this extreme multi-label classification problem (more than
5000 EC numbers in UniProt). Thus, obtaining reliable EC number prediction
results is not straightforward, and the prediction performance is not ideal.
The fourth noteworthy challenge is the usability of existing tools that need
refinement so that the end-user can use them smoothly even with no coding
experience.

In this paper, we take a unified approach to address these challenges. For the
first challenge, we constructed three standard datasets for benchmarking and
evaluation. The datasets contain more than 470,000 distinct labeled protein
sequences from Swiss-Prot. To address the second challenge, we introduced the
cutting-edge ideology from natural language embedding for protein sequence
representation. Firstly, state-of-the-art deep learning methods were evaluated
and adopted for universal protein sequence embedding [21, 22]. Then, we used a
feedback mechanism to choose the most suitable method in response to the
downstream tasks for optimization. To address the third challenge, we proposed
a Dual-core Multiagent Learning Framework (DMLF) for EC number prediction. In
DMLF, we formulate the EC number prediction as a three-step hierarchical
extreme multi-label classification problem. The first step predicts whether a
given protein sequence is an enzyme or not. The second predicts how many
functions the enzyme can perform, i.e., multifunctional enzyme prediction. The
last step predicts the exact EC number for each enzyme function. We use
traditional machine learning methods in the first two steps and a novel deep
learning-based extreme multi-label classifier in the last step, then use a
greedy strategy to integrate these steps to maximize the EC prediction
performance. To address the last challenge, we streamlined the construction
process and open-sourced our codes. Moreover, we published a webserver, so
that anyone can annotate EC numbers smoothly in high-throughput, whether they
have coding experience or not.

##  2 Methodology

This section consists of five subsections. We first formulate the enzyme
function annotation problem in the first subsection, and then describe the
benchmark data construction process in the second subsection. The third
subsection describes our proposed DMLF framework for the benchmark tasks. The
fourth subsection describes the baselines. In the last subsection, we describe
the evaluation metrics.

###  2.1 Problem Formulation

In order to annotate the enzyme function of a new protein sequence, the
initial and basic task is to define whether a given protein is an enzyme.
Since there are numerous multifunctional enzymes, the next task to consider is
to determine whether the enzyme is monofunctional or multifunctional. If it is
multifunctional, the number of functions needs to be classified. After
completing the above two tasks, it is necessary is to assign an EC number to
each function. Based on these considerations, we proposed three basic tasks
for the functional annotation of enzymes, as shown below.

####  2.1.1 Enzyme or Non-enzyme Annotation.

The enzyme or non-enzyme annotation task is formulated as a binary
classification problem:

$f:X\rightarrow\\{0,1\\}$ (1)

where $X=\\{x_{1},x_{2},\cdots,x_{n}\\},n\geq 1$ represents a group of protein
sequences, and $\\{0,1\\}$ is the label indicting whether a given protein is
an enzyme .

####  2.1.2 Multifunctional Enzyme Annotation.

Multifunctional enzyme annotation is formulated as a multi-classification
problem:

$f:X\rightarrow\\{1,2,\cdots,k\\},$ (2)

where $k$ represents the maximum number of EC number for a given protein.

####  2.1.3 Enzyme Commission Number Assignment.

The enzyme commission number assignment task is also formulated as a multi-
classification problem as defined in Eq. 3.

$f:X\rightarrow\\{1.1.1.1,1.1.1.2,1.1.1.3,\cdots\\},$ (3)

###  2.2 Dataset Description

To address the first challenge, we constructed three standard datasets
(Supplemental, SI Appendix Materials and Methods, A. Dataset). Similar to
previous work [21, 26], these datasets are extracted from the Swiss-Prot
database. To simulate real application scenarios as closely as possible, we
did not shuffle data randomly. Instead, after data preprocessing (See
Supplemental, SI Appendix Materials and Methods A. Preprocessing), we
organized data in chronological order. Specifically, we used a snapshot from
Feb 2018 as the training dataset. The training data contains 469,134 distinct
sequences in a total 556,825 records, among which 53.56% are non-enzymes,
while the remaining 47.44% are enzymes. The testing data was extracted from
the June 2020 snapshot and sequences that appeared in the training set were
excluded. The details are listed in Table 1.

Table 1: Description of Benchmarking Datasets | Snapshot | Difference
---|---|---
ITEM |  Feb.2018 |  Jun.2020 |  differ |  added |  deleted
Records |  556,825 |  563,972 |  7147 |  - |  -
Duplicate Removal |  469,134 |  476,006 |  6877 |  8033 |  1156
Non-enzyme |  246,567 |  247,324 |  757 |  4454 |  879
Enzyme |  222,567 |  228,687 |  6120 |  3579 |  277
Distinct EC |  4854 |  5306 |  452 |  644 |  192

$\blacksquare$ Dataset 1: Enzyme and Non-enzyme Dataset

As listed in Table 2, the training set in total has 469,134 records, 222,567
of which are enzymes, and 246,567 are non-enzymes. The testing set contains
7101 records, 3304 of which are enzymes, and the other 3797 are non-enzymes.
To make the data more inclusive, we did not filter any sequence in terms of
length and homology, which is different from previous studies. An enzyme is
labeled as 1 and non-enzyme is labelled as 0. More details about the dataset
can be found in the Supplemental, SI Appendix Materials and Methods, A.
Dataset.

Table 2: Description of Enzyme and Non-enzyme Datasets ITEM |  Training set |  Testing set
---|---|---
Enzyme |  222,567 |  3,304
Non-enzyme |  246,567 |  3,797
Total |  469,134 |  7,101

$\blacksquare$ Dataset 2: Multifunctional Enzyme Dataset

The multifunctional enzyme dataset only contains enzyme data (225,871
records). The number of EC categories ranges from 1 to 8. The details of the
dataset are listed in Table 3.

Table 3: Description of the Multifunctional Enzyme Dataset $\\#$EC | Records | $\\#$EC | Records
---|---|---|---
Trainset | Testset | Trainset | Testset
1 | 210788 | 3052 | 5 | 206 | 6
2 | 9943 | 183 | 6 | 80 | 2
3 | 993 | 53 | 7 | 27 | 1
4 | 525 | 7 | 8 | 5 | 0

$\blacksquare$ Dataset 3: EC number Dataset

Similar to the multifunctional enzyme dataset, the EC number dataset consists
of 225,871 enzyme records, 222,567 of which constitute the training dataset,
and the remaining 3304 are the testing dataset, covering 5111 EC numbers. As
shown in Fig. 1, the test data include 257 newly added EC numbers compared
with the training data, which means that these EC numbers did not appear in
the training process, so predictive methods cannot handle this part of the EC
numbers. Thus, we excluded the sequences with these 257 EC numbers in the
evaluation process. More details about the dataset can be found in
Supplemental, SI Appendix Materials and Methods, A. Dataset.

Figure 1:  Venn diagram of the training and testing datasets

###  2.3 Proposed Framework

To address the second and third challenges: lack of a generic method with high
EC prediction performance and an efficient universal protein sequence
embedding method, we proposed the DMLF approach, composed of an embedding core
and a learning core. These two cores operate relatively independently. The
embedding core is responsible for embedding protein sequences into a machine-
readable matrix. The learning core is responsible for solving specific
downstream biological tasks (e.g., enzyme and non-enzyme prediction,
multifunctional enzyme prediction, and EC number prediction). The overall
scheme of DMLF is illustrated in Fig. 2.

Figure 2:  DMLF is an explicitly designed dual-core driven framework for EC
number prediction. It consists of 2 independent operation units - an embedding
core and a learning core. The embedding core is tasked with converting protein
sequences into features. The learning core is designed to address the specific
biological tasks defined in the problem formulation section. We use different
agents to solve different tasks. Agent 1 was designed to solve the enzyme or
non-enzyme classification task, agent 2 was designed to solve the
multifunctional enzyme prediction task, and agent 3 was designed to solve the
EC number assignment task.

$\blacksquare$ Core 1: Embedding

The objective of this core is to calculate the embedding representations for
protein sequences. For protein sequence encoding/embedding, recent studies
have shown the superior performance of deep learning-based methods compared to
traditional methods [23, 24]. Accordingly, we only compared one-hot encoding
to show the difference between these two kinds of embedding in this study.
Here, we adopted three different embedding methods to calculate the sequence
embedding patterns that adequately represent protein sequences. The first one
is commonly the used one-hot encoding [25]. The second is Unirep [21], an
mLSTM "babbler" deep representation learner for proteins. We used the last
layer for protein representation. The third is the evolutionary scale modeling
embedding method (ESM) [22], a pretrained transformer language model for
protein representation. We used the hidden states from the 1st, 32nd, 33rd
layers as protein embeddings.

$\blacksquare$ Core 2: Learning

The learning core is specialized to perform specific biological tasks using
different agents. In this work, the learning core includes three agents.
Agent-1 is a binary classifier that performs enzyme or non-enzyme prediction.
This classifier was constructed using KNN [26]. Agent-2 is a multi-classifier
that predicts the number of putative functions for a given enzyme. It was
implemented using an integrated sequence aligner, a gradient boost decision
tree, and XGBoost. Agent-3 is also a multi-classifier that performs the EC
number prediction task. As EC number prediction is an extreme multilabel
classification (5852 classes in this benchmark), the performance of
traditional multilabel classification methods such XGBoost, decision tree, and
SVM is abysmal (less than 5% in terms of accuracy). Therefore, we trained a
scalable linear extreme classifier (SLICE)[27] to obtain a more reliable
classification performance in this study. The details of agent implementation
and parameter settings can be found in Supplemental, SI Appendix Materials and
Methods C. Models.

$\blacksquare$ Integration, fine-tuning and output As illustrated in Fig. 2,
the final EC number prediction output is an integrated process. As shown in
Eq. 4, we formulated this integrated process as an optimization problem:

$\centering\mathop{MAX}\limits_{F1}\\{f(ag_{1},ag_{2},ag_{3},sa)\\}\@add@centering$
(4)

where $ag_{1}$, $ag_{2}$, and $ag_{3}$ are the respective prediction results
from Agent-1, Agent-2, and Agent-3, while $sa$ is the predicted result from
multiple sequence alignment. The integration and fine-tuning process aims to
maximize the optimizing objective. In this work, the objective is the
performance of EC number prediction in terms of the F1 score. We used a greedy
strategy to perform this optimization.

###  2.4 Compared Baselines

To evaluate our proposed method comprehensively, we compared our proposed
method with four existing state-of-the-art techniques with ’GOOD’ usability
(see Supplemental, SI RELATED WORK). Four state-of-the-art techniques are:
CatFam, PRIAM (version 2), ECPred, and DeepEC.

###  2.5 Evaluation Metrics

To comprehensively evaluate the proposed method and existing baselines, we use
5 metrics to evaluate binary classification problems and 4 metrics to evaluate
multiple classification problems. For the binary classification task, the
evaluation criteria include ACC(accuracy), PPV (positive predictive value,
precision), NPV(negative predictive value), RC (recall), and F1 value:

$ACC=\frac{TP+TN}{TP+FP+TN+FN+UP+UN}$ (5)

$PPV=\frac{TP}{TP+FP}$ (6)

$NPV=\frac{TN}{TN+FN}$ (7)

$Recall=\frac{TP}{TP+FN+UP}$ (8)

$F1=\frac{2\times PPV\times Recall}{PPV+Recall}$ (9)

where $TP$ is the true positive value that represents the number of samples
correctly identified as positive, $FP$ is the false positive value that
represents the number of samples wrongly identified as positive, $TN$ is the
true negative value that represents the number of samples correctly identified
as negative, $FN$ is false negative value that represents the number of
samples wrongly identified as negative, $UP$ is unclassified positive samples,
and $UN$ is unclassified negative samples.

For multiple classification problems, the evaluation criteria included mACC
(macro-average accuracy), mPR(macro-average precision), mRecall(macro-average
recall), and mF1(macro-average F1 value):

$mACC=\frac{\sum_{i=1}^{n}{ACC_{i}}}{n},{n=1,2,3,\cdots,N}$ (10)

$mPR=\frac{\sum_{i=1}^{n}{PPV_{i}}}{n},{n=1,2,3,\cdots,N}$ (11)

$mRecall=\frac{\sum_{i=1}^{n}{Recall_{i}}}{n},{n=1,2,3,\cdots,N}$ (12)

$mF1=\frac{2\times mPR\times mRecall}{mPR\times mRecall},{n=1,2,3,\cdots,N}$
(13)

where $N$ represents the total number of classes, while $ACC_{i}$, $PPV_{i}$,
and $Recall_{i}$ represent the accuracy, precision, and recall of the $i$-th
class in a one-VS-all mode[28], respectively.

##  3 Results

###  3.1 Embedding Core Performance Evaluation

Figure 3:  Performance comparison of different embedding methods for enzyme or
non-enzyme annotation.

We evaluated five different protein embedding methods, one-hot embedding,
Unirep embedding, and ESM embedding with three different layers (0,32,33) in
our three proposed tasks. We used six machine learning baselines, including
K-nearest neighbor (KNN), logistic regression (LR), XGBoost, decision tree
(DT), random forest (RF), and gradient boosting decision tree (GBDT) to
conduct this evaluation. For embedding, ESM-32 exhibited the best overall
performance among all six baselines regarding all evaluation metrics for
embedding (see Supplemental Tables S12 and S13). As shown in Fig. 3, in task
1, ESM-32 achieved 21.67 and 6.03% improvements over one-hot and Unirep in
terms of accuracy, as well as 27.20 and 7.32% in terms of F1, respectively
(see Supplemental Table S11). This experiment suggests that better embedding
can lead to better learning performance, and deep latent representation can
comprehensively represent the protein sequence. The embedding performance of
ESM-32 was better than that of ESM-33, suggesting that a deeper embedding
layer is not always better. DMLF can automatically choose the best embedding
methods based on the downstream tasks, and ESM-32 exhibited the best
performance in this work.

###  3.2 Task 1: Enzyme or Non-enzyme Prediction

In this work, the workflow of enzyme number assignment is: task 1, determine
whether the given protein sequence is an enzyme $\overrightarrow{}$ task 2, if
the given protein is an enzyme, then predict how many enzyme functions it can
perform; $\overrightarrow{}$ task3, assign an EC number for each enzyme
function. According to this workflow, the first benchmarking task is enzyme or
non-enzyme prediction. In this task, we trained an integrated binary
classification model, which is driven by KNN and sequence alignment. KNN was
implemented using scikit-learn, and the alignment was implemented using
diamond v2.0.11.

Figure 4:  Task 1 Comparison of enzyme or non-enzyme annotation

As shown in Fig. 4, our method can achieve scores of 93.12, 95.25, and 88.99%
in terms of accuracy, precision, and recall, respectively (see Supplemental
Table S8). Compared with other state-of-the-art tools and techniques the
overall accuracy was greatly improved. For example, DeepEC yielded 74.10%,
compared with 93.12% using our algorithm. Many previous methods were designed
to obtain high precision while neglecting accuracy, NPV, and recall. For
example, DeepEC can reach 94.68% precision while recall is only 20.83%.
Methods that only offer high precision are very likely to miss many new
functions. The F1 score might be a better evaluation metric for the EC
assignment of real-world proteins.

###  3.3 Task 2: Multifunctional Enzyme Prediction

The second benchmarking task we addressed is multifunctional enzyme
prediction. The backward prediction engine is agent 2 (see Fig. 2). In this
task, we trained an integrated multiple-classification model driven by
sequence alignment and XGBoost.

Figure 5:  Task 2, Comparison of multifunctional enzyme prediction results.

As shown in Fig. 5, our method was superior to existing baselines (see
Supplemental Table S9). For example, the accuracy and recall of DeepEC were
8.52 and 13.6%, respectively. Moreover, the f1-macro of ECPred, DeepEC,
CatFam, and PRIAM-V2 was less than 6%, lower than a random prediction accuracy
of 10%. Hence, the performance was notably insufficient when dealing with
multifunctional enzyme prediction. The low performance is mainly due to a lack
of multifunctional enzyme data (see Table 3). Although our proposed method is
6.3 times better than random prediction in terms of accuracy, the performance
is still insufficient, so it should be further improved in future work.

###  3.4 Task 3: Enzyme Commission Number Prediction

This task corresponds to agent-3 in DMLF. In order to develop a balanced EC
number prediction algorithm with high accuracy combined with reasonable
precision and recall, we trained an extreme multi-label classification model.

Figure 6:  Task 3 Comparison of EC number prediction results.

Our method achieved 86.91% accuracy with 69% precision and 63.88 recall (Fig.
6), which means that if 100 protein sequences were uploaded for annotation, we
can obtain approximately 87 correct annotations. PRIAM is mainly designed to
include more sequences, so the recall is high (78.48%), while the accuracy
(3.0%) and precision (20.80%) are very low. DeepEC, ECPred, and CATFAM pursue
high precision, so the accuracy is very low (less than 7.5%), which means that
if we upload 100 protein sequences for annotation, we can only obtain 7.5
correct annotations while the remaining 92.5 are wrong. Obviously, our method
shows a clear advantage in terms of EC number assignment.

##  4 Web Server Implementation

To make the workflow accessible for biologists around the world, we built a
web application (https://ecrecer.biodesign.ac.cn/, Fig. 7). End-users can
simply upload sequences to our platform, and then click the submit button to
trigger the prediction workflow. In general, the whole workflow can be
completed in a few seconds. We use Amazon DynamoDB to store job information,
and users can track the previous submission records and corresponding status
information. Once the analysis is finished, the user can view or download the
corresponding results.

For EC assignment workflow, we use Amazon ECR to store Docker images, which
packages a set of bioinformatics software, such as diamond and in-house python
scripts. We built a scalable, elastic, and easily maintainable batch engine
using AWS Batch. This solution took care of dynamically scaling our computer
resources in response to the number of runnable jobs in our job queue.
Finally, we used AWS step functions to coordinate the components of our
applications easily, process message passed from AWS API Gateway, and invoke
the workflows asynchronously. AWS API Gateway was used as the API server to
handle the HTTP requests and route traffic to the correct backends. The static
website was hosted by AWS S3 and sped up using AWS CloudFront.

Figure 7:  The architecture of the web platform.

##  5 Case Study and Discussion

When dealing with the EC assignment problem in a daily production scenario,
ECRECer offers two optional modes for end-users: a prediction mode and a
recommendation mode. In prediction mode, we provide the results with the
highest probability, while in the recommendation mode, we deliver up to 20
possible EC number annotations ranked by their respective likelihood. Here, we
present an up-to-date EC number prediction case to simulate the real-time
challenge by conducting EC number assignments in the prediction mode. We
collected testing protein sequences from June 2020 to November 2021,
encompassing 1968 records. These data were not employed in the development
process of the existing methods or our proposed method, which is in line with
the daily production scenario. In this evaluation, we compared our method with
the state-of-the-art method DeepEC. The comparison results demonstrated the
exceptional performance of our proposed method in terms of accuracy and F1
score. Specifically, our method successfully predicted 1739 records, while
DeepEC correctly predicted 1025 only records (88.36 vs. 52.08%, 60.40 vs.
10.96% in terms of accuracy and F1, respectively, 992 identical prediction
tasks). For the sake of comparison with methods not handling the 7th EC class
(translocases), we reconstructed a reduced sub-dataset comprising 997 enzyme
sequences labeled with EC codes from classes 1 to 6. In this scenario, our
methods still deliver a better performance. ECRECer and DeepEC predicted 807
and 100 correct enzyme records, respectively (96 identical prediction tasks).
This proved the outstanding performance of ECRECer in enzyme annotation.

We present another use case to show that our tools can be used for EC number
completion. In the databases, many enzymes with EC numbers exist in an
uncompleted three-level, two-level, or even one-level state. However, these
proteins with incomplete EC numbers cannot directly be utilized for retrieving
enzymatic reactions. Here, end-users can utilize ECRECer for EC number
completion. For example, in the above case, we mis-predicted 163
monofunctional enzymes, 67 of which had incomplete EC numbers, while ECRECer
completed 38 records to the final fourth level. For instance, for hybrid PKS-
NRPS synthetase Phm1 (UniProt ID: A0A2Z5XAL7, EC: 2.3.1.-), ECRECer predicted
the fourth-level EC number 2.3.1.41. Interestingly, after numerous literature
reviews, we found that this protein has an active site based on PROSITE-
ProRule annotation (PROSITE-ProRule: PRU10022), and it is a beta-ketoacyl
synthases active site (https://prosite.expasy.org/rule/PRU10022). Therefore,
it is very likely that this protein indeed has beta-ketoacyl-acyl-carrier-
protein synthase activity. Another example is iron/alpha-ketoglutarate-
dependent dioxygenase AusU (UniProt ID: A0A0U5GJ41). This protein has a two-
level EC number in the database (1.14.-.-). When we used ECRECer for EC
annotation, it assigned this protein with the fourth-level EC 1.14.11.38. This
protein was recently integrated into UniProtKB/Swiss-Prot (September 29,
2021). After blasting it against the UniProtKB database, we found that the top
5 reviewed proteins with the highest identities include three verruculogen
synthase (Fig. 8a). We took protein Q4WAW9 as an example, and found that both
genes belong to exactly the same protein families with the same domains (Fig.
8b). To further validate the results, we compared the structure of A0A0U5GJ41
(alphfold2 predicted) and Q4WAW9 (alphfold2 predicted and crystal structure).
The results showed that these two proteins have a highly similar structure
(see Supplemental Figs. S3-S6) with small RMSD (1.104). Therefore, the protein
could be potentially annotated as EC 1.14.11.38 as well.

Figure 8:  a) Blast search of the protein sequence against the UniProtKB
database; b) Annotation of protein families and domains; c)Comparison of
structural similarity.

In addition to EC number assignment, another advantage of ECRECer is the
recommendation of EC numbers, which makes our tool unique. The recommendation
is particularly helpful for the discovery of multifunctional enzymes. For
example, in the case of annotating the intradiol ring-cleavage dioxygenase
PrcA (UniProt ID: A2R1P9, EC: 1.13.11.3), the prediction mode indicates
EC:1.13.11.1 as the most likely EC number, while the recommendation mode
presents a list of recommendations [1.13.11.1, 1.13.11.3, 3.5.2.17, … ]. The
correct EC number is the second in this list. Actually, the substrates of the
two reactions are very similar, and 3,4-dihydroxybenzoate only has an
additional carboxyl group compared with catechol. Hence, PrcA might be a
multifunctional enzyme. Another example is laccase (UniProt ID: A0A7T1FRB0,
EC: 1.10.3.2), as its functional annotation is relatively sparce in the
reviewed Swiss-Prot data. In this case, DeepEC, sequence alignment, and our
prediction mode give out the prediction of a non-enzyme, but when we used the
recommendation mode, we obtained a recommendation list [non-enzyme, 1.10.3.2,
1.-.-.-, 6.1.1.3, …], and the correct annotation as offered as the second
recommendation.

To demonstrate the inclusiveness and predictive ability of our proposed
method, we conducted EC number prediction on an unreviewed protein family.
Corynebacterium glutamicum, the famous industrial workhorse for amino acid
production with a current output of over 6 million tons per year (Lee et al.,
2016), is increasingly being adopted as a promising chassis for the
biosynthesis of other compounds. However, unlike E. coli (1652 protein
sequences with EC numbers out of 4322 proteins, 38.2%), the protein sequences
of Corynebacterium glutamicum were not well annotated. Out of 3305 protein
sequences, only 537 were reviewed and included in the Swiss-Prot database (357
proteins have assigned EC numbers). We used the other 2768 protein sequences
to compare our tool with DeepEC. Our approach was able to assign 1056 proteins
with EC numbers s, while DeepEC only assigned 157 EC numbers (123 same EC
numbers between DeepEC and ECRECer). Although there is no gold standard to
decide which prediction is correct, we believe our algorithm should provide a
more reasonable prediction as the proportion of protein sequences with EC
numbers is similar to that of E. coli (42% vs. 15.5% in the case of DeepEC).
The newly predicted EC numbers for the protein sequences are crucial for
further analysis, such as retrieving metabolic reactions for genome-scale
modeling.

##  6 Conclusion

In this work, we proposed a novel dual-core multiagent learning framework to
complete three benchmarking tasks: 1) enzyme or non-enzyme annotation; 2)
predicting the number of catalytic functions of a single multifunctional
enzyme; and 3) EC number prediction. The method developed in this work has two
calculation cores, an embedding core and a learning core. The embedding core
is responsible for selecting the best available embedding method among one-
hot, Unirep, and ESM to calculate sequence embeddings. The learning core is
responsible for completing the specific benchmarking tasks using the best
calculated protein sequence embedding as input.

We were guided by two principles in the design of our methods. The first
principle is high usability (both can be accessed via the world-wide-web and
provide standalone suit for high throughput prediction) with relatively
balanced prediction performance (which can achieve the best accuracy with
reasonable precision and recall). The second principle is providing
comprehensive evaluation metrics with accessible reproduction datasets and
source codes. To implement the first principle, we proposed DMLF. To implement
the second principle, we provided a web server, standalone packages and opened
all the source codes, including data preprocessing, dataset buildup, model
training, and model testing/evaluation.

Experiments on real-world datasets and comprehensive comparisons with existing
state-of-the-art methods demonstrated that our tool is highly competitive, has
the best performance with high usability, and meets the proposed objectives.
Although our tool exhibited the best performance, it still has much space for
improvement. For example, the performance of multifunctional enzyme annotation
is relatively low, while the accuracy and recall of EC number annotation is
less than 90%. Our feature work will focus on improving the prediction
precision.  Key Points • A multiagent dual-core learning framework is proposed
to predict Enzyme Commission (EC) Numbers by using protein sequence data. • A
protein language model and an extreme multi-label classifier are adopted to
reduce the heavy head-crafted feature engineering and elevate the prediction
performance. • The proposed framework remarkably outperforms the existing
state-of-the-the-art method in terms of accuracy and F1 score by 70% and 20%,
respectively. • An online service and an offline bundle are provided for end-
users to annotate EC numbers in high-throughput easily and efficiently.

##  7 Supplementary Data

Supplementary data are available onle at
https://github.com/kingstdio/ECRECer/blob/main/document/supplementary.pdf

##  8 Data availability

The data underlying this article are available in the article and in its
online Supplementary Material. The code of ECRECer,the training data,and the
prediction results areavailable at urlhttps://ecrecer.biodesign.ac.cn/.

##  9 Author contributions statement

Z.S. and X.L. designed and implemented the model, conducted the experiments,
analyzed the results and wrote the manuscript. H.M. and Q.Y. reviewed the
manuscript. R.W. and H.L. designed the website.

##  10 Funding

This work was supported by the National Key Research and Development Program
of China [2018YFA0900300, 2020YFA0908300], the International Partnership
Program of Chinese Academy of Sciences [153D31KYSB20170121], the Youth
Innovation Promotion Association CAS, and Tianjin Synthetic Biotechnology
Innovation Capacity Improvement Project [TSBICIP-PTJS-001, TSBICIP-CXRC-018].

## References

  * [1] The UniProt Consortium.  Uniprot: the universal protein knowledgebase in 2021.  Nucleic Acids Research, 49(D1):D480–D489, 2021.
  * [2] Jae Yong Ryu, Hyun Uk Kim, and Sang Yup Lee.  Deep learning enables high-quality and high-throughput prediction of enzyme commission numbers.  Proceedings of the National Academy of Sciences, 116(28):13996–14001, 2019.
  * [3] Nicholas Furnham, John S Garavelli, Rolf Apweiler, and Janet M Thornton.  Missing in action: enzyme functional annotations in biological databases.  Nature chemical biology, 5(8):521–525, 2009.
  * [4] Wikimedia.  Enzyme commission number.  [EB/OL], 2021.  https://en.wikipedia.org/wiki/Enzyme_Commission_number Accessed November 29, 2021.
  * [5] Jui-Hung Hung and Zhiping Weng.  Sequence alignment and homology search with blast and clustalw.  Cold Spring Harbor Protocols, 2016(11):pdb–prot093088, 2016.
  * [6] Chenggang Yu, Nela Zavaljevski, Valmik Desai, and Jaques Reifman.  Genome-wide enzyme annotation with precision control: Catalytic families (catfam) databases.  Proteins: Structure, Function, and Bioinformatics, 74(2):449–460, 2009.
  * [7] Clotilde Claudel-Renard, Claude Chevalet, Thomas Faraut, and Daniel Kahn.  Enzyme-specific profiles for genome annotation: Priam.  Nucleic acids research, 31(22):6633–6639, 2003.
  * [8] Nirvana Nursimulu, Leon L Xu, James D Wasmuth, Ivan Krukov, and John Parkinson.  Improved enzyme annotation with ec-specific cutoffs using detect v2.  Bioinformatics, 34(19):3393–3395, 2018.
  * [9] Chengxin Zhang, Peter L Freddolino, and Yang Zhang.  Cofactor: improved protein function prediction by combining structure, sequence and protein–protein interaction information.  Nucleic acids research, 45(W1):W291–W299, 2017.
  * [10] Ying Hong Li, Jing Yu Xu, Lin Tao, Xiao Feng Li, Shuang Li, Xian Zeng, Shang Ying Chen, Peng Zhang, Chu Qin, Cheng Zhang, et al.  Svm-prot 2016: a web-server for machine learning prediction of protein functional families from sequence irrespective of similarity.  PloS one, 11(8):e0155290, 2016.
  * [11] Alperen Dalkiran, Ahmet Sureyya Rifaioglu, Maria Jesus Martin, Rengul Cetin-Atalay, Volkan Atalay, and Tunca Doğan.  Ecpred: a tool for the prediction of the enzymatic functions of protein sequences based on the ec nomenclature.  BMC bioinformatics, 19(1):1–13, 2018.
  * [12] Adrian K Arakaki, Ying Huang, and Jeffrey Skolnick.  Eficaz 2: enzyme function inference by a combined approach enhanced by machine learning.  BMC bioinformatics, 10(1):1–15, 2009.
  * [13] Taofeek D Akinosho, Lukumon O Oyedele, Muhammad Bilal, Anuoluwapo O Ajayi, Manuel Davila Delgado, Olugbenga O Akinade, and Ashraf A Ahmed.  Deep learning in the construction industry: A review of present status and future innovations.  Journal of Building Engineering, page 101827, 2020.
  * [14] Haoyang Li, Shuye Tian, Yu Li, Qiming Fang, Renbo Tan, Yijie Pan, Chao Huang, Ying Xu, and Xin Gao.  Modern deep learning in bioinformatics.  Journal of molecular cell biology, 12(11):823–827, 2020.
  * [15] Yuxi Li, Yue Zuo, Houbing Song, and Zhihan Lv.  Deep learning in security of internet of things.  IEEE Internet of Things Journal, 2021.
  * [16] Zhenkun Shi, Sen Wang, Lin Yue, Lixin Pang, Xianglin Zuo, Wanli Zuo, and Xue Li.  Deep dynamic imputation of clinical time series for mortality prediction.  Information Sciences, 579:607–622, 2021.
  * [17] Hong-Bin Shen and Kuo-Chen Chou.  Ezypred: a top–down approach for predicting enzyme functional classes and subclasses.  Biochemical and biophysical research communications, 364(1):53–59, 2007.
  * [18] Yu Li, Sheng Wang, Ramzan Umarov, Bingqing Xie, Ming Fan, Lihua Li, and Xin Gao.  Deepre: sequence-based enzyme ec number prediction by deep learning.  Bioinformatics, 34(5):760–769, 2018.
  * [19] Ji Yong An, Yong Zhou, Yu Jun Zhao, and Zi Ji Yan.  An Efficient Feature Extraction Technique Based on Local Coding PSSM and Multifeatures Fusion for Predicting Protein-Protein Interactions.  Evolutionary Bioinformatics, 15, 2019.
  * [20] Kevin K. Yang, Zachary Wu, Claire N. Bedbrook, and Frances H. Arnold.  Learned protein embeddings for machine learning.  Bioinformatics, 34(15):2642–2648, 2018.
  * [21] Ethan C Alley, Grigory Khimulya, Surojit Biswas, Mohammed AlQuraishi, and George M Church.  Unified rational protein engineering with sequence-based deep representation learning.  Nature methods, 16(12):1315–1322, 2019.
  * [22] Roshan Rao, Joshua Meier, Tom Sercu, Sergey Ovchinnikov, and Alexander Rives.  Transformer protein language models are unsupervised structure learners.  In International Conference on Learning Representations, 2020.
  * [23] Marco Anteghini, Vitor Martins dos Santos, and Edoardo Saccenti.  In-pero: Exploiting deep learning embeddings of protein sequences to predict the localisation of peroxisomal proteins.  International Journal of Molecular Sciences, 22(12):6409, 2021.
  * [24] Hannah-Marie Martiny, Jose Juan Almagro Armenteros, Alexander Rosenberg Johansen, Jesper Salomon, and Henrik Nielsen.  Deep protein representations enable recombinant protein expression prediction.  bioRxiv, 2021.
  * [25] Hesham ElAbd, Yana Bromberg, Adrienne Hoarfrost, Tobias Lenz, Andre Franke, and Mareike Wendorff.  Amino acid encoding for deep learning applications.  BMC bioinformatics, 21(1):1–14, 2020.
  * [26] Shichao Zhang, Xuelong Li, Ming Zong, Xiaofeng Zhu, and Debo Cheng.  Learning k for knn classification.  ACM Transactions on Intelligent Systems and Technology (TIST), 8(3):1–19, 2017.
  * [27] Himanshu Jain, Venkatesh Balasubramanian, Bhanu Chunduri, and Manik Varma.  Slice: Scalable linear extreme classifiers trained on 100 million labels for related searches.  In Proceedings of the Twelfth ACM International Conference on Web Search and Data Mining, pages 528–536, 2019.
  * [28] Ryan Rifkin and Aldebaro Klautau.  In defense of one-vs-all classification.  The Journal of Machine Learning Research, 5:101–141, 2004.


# The average distance problem with perimeter-to-area ratio
penalization111This paper will appear in SIAM Journal on Mathematical
Analysis.

Qiang Du Department of Applied Physics and Applied Mathematics, and Data
Science Institute, Columbia University, 500 W. 120th St., New York, NY 10027,
USA. Email<EMAIL_ADDRESS>Xin Yang Lu Department of Mathematical
Sciences, Lakehead University, 955 Oliver Rd., Thunder Bay, ON P7B 5E1, Canada
AND Department of Mathematics and Statistics, McGill University, 805
Sherbrooke St. W., Montreal, QC H3A 0B9, Canada. Email<EMAIL_ADDRESS>Chong Wang Department of Mathematics, Washington and Lee University, 204 W
Washington St., Lexington, VA 24450, USA. Email<EMAIL_ADDRESS>
###### Abstract

In this paper we consider the functional

$E_{p,\lambda}(\Omega):=\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}.$

Here $p\geq 1$, $\lambda>0$ are given parameters, the unknown $\Omega$ varies
among compact, convex, Hausdorff two-dimensional sets of $\mathbb{R}^{2}$,
$\partial\Omega$ denotes the boundary of $\Omega$, and
$\operatorname{dist}(x,\partial\Omega):=\inf_{y\in\partial\Omega}|x-y|$. The
integral term
$\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
quantifies the “easiness” for points in $\Omega$ to reach the boundary, while
$\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}$ is the
perimeter-to-area ratio. The main aim is to prove existence and
$C^{1,1}$-regularity of minimizers of $E_{p,\lambda}$.

Keywords. perimeter-to-area ratio, regularity

Classification. 49Q20, 49K10, 35B65

##  1 Introduction

The perimeter-to-area ratio (in 2D), or surface area-to-volume ratio (in 3D),
plays a crucial role in many processes. In biology, for instance, the size of
prokaryote cells is limited by the efficiency of diffusion processes,
fundamental to transport nutrients across the cell, which is strongly
correlated with the surface area-to-volume ratio. A larger surface area-to-
volume ratio also gives prokaryote cells a high metabolic rate, fast growth,
and short lifespan compared to eukaryote cells (see for instance [10]).

In chemistry, higher surface area-to-volume ratio increases the typical speed
of chemical reactions. This phenomenon can be observed in many instances ,
sometimes quite dramatically, such as dust explosions, when dust particles of
seemingly non-flammable materials (e.g., aluminum, sugar, flour, etc.) can be
ignited due to their very large surface area-to-volume ratio ([14, 12]).

In this paper we will focus on the 2D case. In the above examples, there are
essentially two often competing quantities: one is the “easiness” to access
the boundary, and the other is the perimeter-to-area ratio.

A very thin, rod-like, rectangular body would have very good access to
boundary (desirable), but large perimeter-to-area ratio. A disk would have the
lowest perimeter-to-area ratio (desirable) among shapes of the same total
area, but access to boundary would be limited. It is also possible to have
both large perimeter-to-area ratio and limited access to boundary.

Until now, we have discussed the “easiness” of accessing the boundary only at
a qualitative level. In order to quantify it, we introduce the “average
distance” term

$F_{p}(\Omega):=\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x,$

where $\operatorname{dist}(x,\partial\Omega):=\inf_{y\in\partial\Omega}|x-y|$;
$p\geq 1$ is a given parameter; and $|\cdot|$ denotes the Euclidean distance.

Consider the energy functional

$\displaystyle E_{p,\lambda}(\Omega)=\int_{\Omega}\text{
dist}^{p}(x,\partial\Omega)\
dx+\lambda\frac{{\mathcal{H}}^{1}(\partial\Omega)}{{\mathcal{H}}^{2}(\Omega)},$
(1.1)

where $p\geq 1,\lambda>0$ are given parameters. Define the admissible set

$\displaystyle{\mathcal{A}}:=\\{\Omega:\Omega\subset\mathbb{R}^{2}\text{ is
compact, convex and Hausdorff two-dimensional}\\}.$

The term $\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}$ is
the perimeter-to-area ratio. Note that neither the perimeter
$\mathcal{H}^{1}(\partial\Omega)$, nor the area $\mathcal{H}^{2}(\Omega)$, is
penalized, only their ratio is. This makes compactness results quite
challenging to prove, and several estimates (in Section 2) will be required.
Another issue is that it is not very clear if the average-distance term is
just a lower order perturbation of
$\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}$. The role of
convexity is to ensure crucial compactness estimates (Lemmas 2.2 and 3.5).
Note that $E_{p,\lambda}$ is invariant under rigid movements. Further details
about the space of convex sets, and its topology, will be discussed in Section
2. The main result of this paper is:

######  Theorem 1.1.

Given $p\geq 1,\lambda>0$, the following assertions hold:

  * (1)

$E_{p,\lambda}$ admits a minimizer in ${\mathcal{A}}$.

  * (2)

All minimizers are compact, convex, $C^{1,1}$-regular sets, with Hausdorff
dimension equal to $2$.

  * (3)

The perimeter-to-area ratio of any minimizer $\Omega$ satisfies

$\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}=\frac{p+2}{\lambda(p+3)}\min_{\mathcal{A}}E_{p,\lambda}.$

Here, and for future reference, the expression “$\Omega$ is $C^{k}$-regular”
means that its boundary $\partial\Omega$ is $C^{k}$-regular, i.e.,
$\partial\Omega$ admits a $C^{k}$-regular parameterization.

Note that the functional $F_{p}$ is formally similar to the average-distance
functional

$\Sigma\mapsto\int_{\Gamma}\operatorname{dist}^{p}(x,\Sigma)\,{\operatorname{d}}\mu,$

where $\Gamma$ is a given domain, $\mu$ a given measure on $\Gamma$, and
$\Sigma$ varies among compact, path-wise connected sets with Hausdorff
dimension equal to 1. The average-distance functional has been widely studied,
and used in several modeling problems. For a (non exhaustive) list of
references, we cite the papers (and books) by Buttazzo and collaborators [2,
3, 4, 8, 9, 6, 7]. Also related are the papers by Paolini and Stepanov [22],
Santambrogio and Tilli [23], Tilli [26], Lemenant and Mainini [19], Slepčev
[25], and the review paper by Lemenant [18]. Similar variational problems
entailing a competition between classical perimeter and nonlocal repulsive
interaction were studied by Muratov and Knüpfer [21], Goldman, Novaga and
Ruffini [16], and Goldman, Novaga and Röger [15]. Figalli, Fusco, Maggi,
Millot, and Morrini studied a competition between a nonlocal $s$-perimeter and
a nonlocal repulsive interaction term [13].

The rest of the paper is structured as follows: section 2 is dedicated to
proving some auxiliary estimates on the area ((2.2) and Corollary 2.1) and
perimeter (Lemma 2.2) of elements of minimizing sequences. Existence of
minimizers will be shown in section 3, while $C^{1,1}$ regularity will be
proven in section 4. Finally, we explore several future directions to further
our understanding of the penalized average distance problem.

##  2 Preliminary estimates

In this section we collect some preliminary estimates that will be used later.
First, we remark that given $p\geq 1$ and $\lambda>0$, for any
$\Omega\in{\mathcal{A}}$ it holds

$\mathcal{H}^{2}(\Omega)\geq\frac{4\pi\lambda^{2}}{E_{p,\lambda}(\Omega)^{2}}.$
(2.2)

Indeed, consider an arbitrary $\Omega\in{\mathcal{A}}$. By the isoperimetric
inequality, among all convex sets with area $\mathcal{H}^{2}(\Omega)$, the
perimeter-to-area ratio is minimum for a disk, where it attains the value
$2\sqrt{\pi}/\sqrt{\mathcal{H}^{2}(\Omega)}$. Hence

$\frac{2\lambda\sqrt{\pi}}{\sqrt{\mathcal{H}^{2}(\Omega)}}\leq\lambda\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\leq
E_{p,\lambda}(\Omega),$

and (2.2) is proven.

######  Corollary 2.1.

Given $p\geq 1$, $\lambda>0$, any minimizing sequence
$\Omega_{n}\subseteq{\mathcal{A}}$ satisfies

$\displaystyle\mathcal{H}^{2}(\Omega_{n})$ $\displaystyle\geq
4\pi\lambda^{2}\bigg{(}\frac{2\pi}{p^{2}+3p+2}+2\lambda+1\bigg{)}^{-2}=:C_{1},$
(2.3)
$\displaystyle\frac{\mathcal{H}^{1}(\partial\Omega_{n})}{\mathcal{H}^{2}(\Omega_{n})}$
$\displaystyle\leq\frac{1}{\lambda}\bigg{(}\frac{2\pi}{p^{2}+3p+2}+2\lambda+1\bigg{)}=:C_{2},$
(2.4)

for any sufficiently large $n$.

###### Proof.

First we prove $\inf_{\mathcal{A}}E_{p,\lambda}<+\infty$. Let
$B_{1}\in{\mathcal{A}}$ be a disk of radius $1$. Direct computation gives

$\displaystyle\inf_{{\mathcal{A}}}E_{p,\lambda}\leq E_{p,\lambda}(B_{1})$
$\displaystyle=\int_{B_{1}}\operatorname{dist}^{p}(x,\partial
B_{1})\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial
B_{1})}{\mathcal{H}^{2}(B_{1})}$
$\displaystyle=2\pi\int_{0}^{1}(1-r)^{p}r\,{\operatorname{d}}r+2\lambda=\frac{2\pi}{p^{2}+3p+2}+2\lambda<+\infty.$
(2.5)

Thus, given a minimizing sequence $\Omega_{n}\subseteq{\mathcal{A}}$, there
exists $N$ such that for any $n\geq N$ it holds

$E_{p,\lambda}(\Omega_{n})\leq\frac{2\pi}{p^{2}+3p+2}+2\lambda+1,$ (2.6)

and (2.2) gives

$\mathcal{H}^{2}(\Omega_{n})\geq
4\pi\lambda^{2}\bigg{(}\frac{2\pi}{p^{2}+3p+2}+2\lambda+1\bigg{)}^{-2},$

for any $n\geq N$, hence (2.3). To prove (2.4), note that (2.6) forces

$\frac{2\pi}{p^{2}+3p+2}+2\lambda+1\geq
E_{p,\lambda}(\Omega_{n})\geq\lambda\frac{\mathcal{H}^{1}(\partial\Omega_{n})}{\mathcal{H}^{2}(\Omega_{n})},$

concluding the proof. ∎

######  Lemma 2.2.

Given $p\geq 1$ and $\lambda>0$, for any minimizing sequence
$\Omega_{n}\subseteq{\mathcal{A}}$, it holds, for all sufficiently large $n$,

$\mathcal{H}^{1}(\partial\Omega_{n})\leq C_{3}=C_{3}(p,\lambda)$ (2.7)

with $C_{3}$ being some computable (but uninfluential) constant.

###### Proof.

We first claim that for any $\Omega\in{\mathcal{A}}$ it holds

$\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x\geq
C\frac{\mathcal{H}^{2}(\Omega)^{p+1}}{\mathcal{H}^{1}(\partial\Omega)^{p}},\qquad
C=3^{-p}2^{-p-4}.$ (2.8)

Consider an arbitrary $\Omega\in{\mathcal{A}}$. Let $A,B\in\partial\Omega$ be
two points realizing $D:=|A-B|=\operatorname{{\operatorname{diam}}}\Omega$.
Let $\Sigma_{i}$, $i=1,2$ be the lines (see Figure 1) orthogonal to the line
segment between $A$ and $B$ (which we denote by $\llbracket A,B\rrbracket$).
Since $\Omega$ is convex, and
$|A-B|=\operatorname{{\operatorname{diam}}}\Omega$, $\Omega$ is entirely
contained in the region between $\Sigma_{1}$ and $\Sigma_{2}$. Then let
$P_{i}$, $i=1,2$ be the points on $\partial\Omega$ such that the triangles
$\triangle AP_{i}B$ have maximal areas. As $\Omega$ is convex, we have

$\mathcal{H}^{2}(\Omega)\leq D(h_{1}+h_{2}),\qquad
h_{i}:=\operatorname{dist}(P_{i},\llbracket A,B\rrbracket).$

On the other hand, $\mathcal{H}^{2}(\triangle AP_{i}B)=Dh_{i}/2$, hence

$\frac{\mathcal{H}^{2}(\triangle AP_{1}B\cup\triangle
AP_{2}B)}{\mathcal{H}^{2}(\Omega)}\geq\frac{1}{2}.$

Figure 1:  A schematic representation of the construction. The points $P_{1}$
and (resp. $P_{2}$) are the points on $\partial\Omega$ above (resp. below) the
segment $\llbracket A,B\rrbracket$ furthest away from $\llbracket
A,B\rrbracket$.

Now we do the following construction: let $O_{i}$ (resp. $r_{i}$) be the
incenter (resp. inradius) of $\triangle AP_{i}B$, $i=1,2$. Denote by
$\tilde{A}$ (resp. $\tilde{P}_{i}$, $\tilde{B}$) the midpoints of the line
segments between $O_{i}$ and $A$ (resp. $P_{i}$, $B$) – see Figure 2.

Figure 2: A schematic representation of the construction. The points
$\tilde{A}$, $\tilde{P}_{i}$, $\tilde{B}$ are the midpoints of the segments
$\llbracket O_{i},A\rrbracket$, $\llbracket O_{i},P_{i}\rrbracket$,
$\llbracket O_{i},B\rrbracket$ respectively. The red dotted circle is the
incircle of the triangle $\triangle AP_{i}B$.

Clearly, $\triangle\tilde{A}\tilde{P}_{i}\tilde{B}$ is a rescaled copy of
$\triangle AP_{i}B$, with area $\mathcal{H}^{2}(\triangle AP_{i}B)/4$. As

$\text{inradius}=\frac{2\text{Area}}{\text{Perimeter}},$

we can estimate $r_{i}$ as follows:

$\displaystyle r_{i}$
$\displaystyle=\frac{Dh_{i}}{D+|A-P_{i}|+|B-P_{i}|}\geq\frac{Dh_{i}}{3D}=\frac{h_{i}}{3}$
(2.9)

since by definition we have
$D=\operatorname{{\operatorname{diam}}}\Omega\geq|A-P_{i}|,|B-P_{i}|$. Then,
noting that

$\operatorname{dist}(x,\partial\Omega)\geq\operatorname{dist}(x,\partial
AP_{i}B)\geq\frac{1}{2}\operatorname{dist}(O_{i},\partial
AP_{i}B)\geq\frac{1}{2}r_{i}$

for all $x\in\triangle\tilde{A}\tilde{P}_{i}\tilde{B},\;i=1,2,$ we have

$\displaystyle\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
$\displaystyle\geq\sum_{i=1}^{2}\int_{\triangle\tilde{A}\tilde{P}_{i}\tilde{B}}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
$\displaystyle\geq\sum_{i=1}^{2}2^{-p}r_{i}^{p}\mathcal{H}^{2}(\triangle\tilde{A}\tilde{P}_{i}\tilde{B})=\sum_{i=1}^{2}2^{-p-2}r_{i}^{p}\mathcal{H}^{2}(\triangle
AP_{i}B)$ $\displaystyle\geq\sum_{i=1}^{2}3^{-p}2^{-p-3}h_{i}^{p+1}D\geq
3^{-p}2^{-p-3}D\cdot\max_{i=1,2}h_{i}^{p+1}.$ (2.10)

Recalling that $\mathcal{H}^{2}(\Omega)\leq D(h_{1}+h_{2})$,
$\mathcal{H}^{1}(\partial\Omega)\geq 2D$, we get

$\frac{\mathcal{H}^{2}(\Omega)^{p+1}}{\mathcal{H}^{1}(\partial\Omega)^{p}}\leq\frac{D^{p+1}(h_{1}+h_{2})^{p+1}}{(2D)^{p}}=2^{-p}D(h_{1}+h_{2})^{p+1}\leq
2D\cdot\max_{i=1,2}h_{i}^{p+1}.$

Hence (2.10) gives

$\displaystyle\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
$\displaystyle\geq
C\frac{\mathcal{H}^{2}(\Omega)^{p+1}}{\mathcal{H}^{1}(\partial\Omega)^{p}},\qquad
C=3^{-p}2^{-p-4},$

and (2.8) is proven. From (2.4) we know that

$\frac{\mathcal{H}^{1}(\partial\Omega_{n})}{\mathcal{H}^{2}(\Omega_{n})}\leq
C_{2}\Longrightarrow\frac{\mathcal{H}^{2}(\Omega_{n})^{p+1}}{\mathcal{H}^{1}(\partial\Omega_{n})^{p+1}}\geq
C_{2}^{-p-1},$

so the above inequality gives

$\displaystyle\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x$
$\displaystyle\geq
C\frac{\mathcal{H}^{2}(\Omega_{n})^{p+1}}{\mathcal{H}^{1}(\partial\Omega_{n})^{p+1}}\mathcal{H}^{1}(\partial\Omega_{n})\geq
CC_{2}^{-p-1}\mathcal{H}^{1}(\partial\Omega_{n}).$

Now, any minimizing sequence $\\{\Omega_{n}\\}$ is such that, for all
sufficiently large $n$,

$E_{p,\lambda}(\Omega_{n})\leq\inf E_{p,\lambda}+1,$

thus

$\displaystyle\inf E_{p,\lambda}+1$ $\displaystyle\geq
E_{p,\lambda}(\Omega_{n})\geq\int_{\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x\geq
CC_{2}^{-p-1}\mathcal{H}^{1}(\partial\Omega_{n}),$

(2.5) shows that $\inf E_{p,\lambda}<+\infty$, completing the proof.

∎

######  Remark 2.3.

We note that it is an interesting geometric question by itself to study what
the optimal constant $C$ for the inequality (2.8). Furthermore, one may ask if
the form of the inequality is optimal. That is, one may ask, given
$\mathcal{H}^{2}(\Omega)$ and $\mathcal{H}^{1}(\partial\Omega)$, what is the
minimum of
$\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$,
which is a constrained optimization problem related to the one considered in
this work.

##  3 Existence

In this section we will prove that the $E_{p,\lambda}$ admits a minimizer in
$\mathcal{A}$. As our arguments rely on a lower semicontinuity result, namely
Lemma 3.4 below, we need first to introduce a metric on $\mathcal{A}$.

For any $\Omega_{1},\Omega_{2}\in{\mathcal{A}}$, define

$\displaystyle
d(\Omega_{1},\Omega_{2}):={\mathcal{H}}^{2}(\Omega_{1}\triangle\Omega_{2}),$
(3.11)

where $\triangle$ denotes the symmetric difference. Set

$\bar{\mathcal{A}}:=\text{completion of }{\mathcal{A}}\text{ with respect to
}d.$

Before we can proceed, we need to characterize the elements of
$\bar{\mathcal{A}}\setminus\mathcal{A}$: we cannot exclude a priori that an
element $\Omega\in\bar{\mathcal{A}}$ can be quite irregular:

  1. 1.

$\Omega\in\bar{\mathcal{A}}$ needs not to be closed: indeed it is very
possible for a sequence of compact sets to converge to an open set in the
metric $d$. For instance, let $\Omega_{n}$ be the closed ball of radius
$1-1/n$ centered around the origin, then it converges to the open ball,
centered around the origin, of radius 1.

  2. 2.

As we do not have any a priori bounds on the diameter of elements of
$\mathcal{A}$, a set $\Omega\in\bar{\mathcal{A}}$ needs not to be bounded.

  3. 3.

The distance $d$ is insensitive to perturbations on
$\mathcal{H}^{2}$-negligible sets. Therefore, we cannot exclude that
$\bar{\mathcal{A}}$ might contain compact convex sets up to
$\mathcal{H}^{2}$-negligible sets. Thus whether a generic element in
$\bar{\mathcal{A}}$ is convex or not is unclear.

In view of the above mentioned issues, we cannot assume neither compactness,
nor convexity, for elements of $\bar{\mathcal{A}}$. Our goal is to show (see
Lemma 3.3 below) that minimizing sequences must converge to some element in
$\mathcal{A}$.

The next result, from [25], will be crucial for our convergence arguments.

######  Lemma 3.1.

Consider a sequence of constant speed parameterized curves
$\gamma_{n}:[0,1]\longrightarrow K$, where $K\subseteq\mathbb{R}^{d}$ is some
compact set. Assume moreover that

$\displaystyle\sup_{n}L(\gamma_{n})<+\infty,\qquad\sup_{n}\|\gamma_{n}\|_{BV([0,1];\mathbb{R}^{d})}<+\infty,$
(3.12)

where $\|\cdot\|_{BV([0,1];\mathbb{R}^{d})}$ denotes the bounded variation
norm. Then there exists a curve $\gamma:[0,1]\longrightarrow K$ such that:

  1. 1.

$\gamma_{n}\to\gamma$ in $C^{\alpha}([0,1];\mathbb{R}^{d})$ for all
$\alpha\in[0,1)$,

  2. 2.

$\gamma_{n}^{\prime}\to\gamma^{\prime}$ in $L^{p}(0,1;\mathbb{R}^{d})$ for all
$p<+\infty$,

  3. 3.

$\gamma_{n}^{\prime\prime}\;{\overset{*}{\rightharpoonup}}\;\gamma^{\prime\prime}$
weakly as measures.

######  Remark 3.2.

We remark that this convergence result is quite strong: consider a sequence
$\\{\Omega_{n}\\}\subseteq\mathcal{A}$ and let $\gamma_{n}$ be constant speed
parameterizations of $\partial\Omega_{n}$. Note that $\gamma_{n}$ are all
closed curves. Assume that we are under the hypotheses of Lemma 3.1, hence
there exists $\gamma:[0,1]\longrightarrow K$ such that $\gamma_{n}\to\gamma$
in $C^{\alpha}([0,1];\mathbb{R}^{d})$ for all $\alpha\in[0,1)$. In particular,
we can define $\Omega$ to be the bounded region delimited by the graph of
$\gamma$, and we have the uniform convergence of the boundaries, which in turn
gives $d_{\mathcal{H}}(\partial\Omega_{n},\partial\Omega)\to 0$. Here
$d_{\mathcal{H}}$ denotes the Hausdorff distance

$d_{\mathcal{H}}(X,Y):=\max\bigg{\\{}\sup_{x\in
X}\operatorname{dist}(x,Y),\sup_{y\in Y}\operatorname{dist}(y,X)\bigg{\\}}.$

Such strong convergence also implies that the characteristic functions
$\chi_{\Omega_{n}}$ converge to $\chi_{\Omega}$ in $L^{p}$, $p\in[1,+\infty)$,
since

$\displaystyle\|\chi_{\Omega_{n}}-\chi_{\Omega}\|_{L^{p}(\mathbb{R}^{d})}^{p}\leq\mathcal{H}^{2}(\Omega_{n}\triangle\Omega)\leq\max\Big{\\{}\sup_{n}L(\gamma_{n}),L(\gamma)\Big{\\}}\cdot
d_{\mathcal{H}}(\partial\Omega_{n},\partial\Omega)\to 0.$

######  Lemma 3.3.

Consider a minimizing sequence $\Omega_{n}\subseteq\mathcal{A}$, then there
exists $\Omega\in\mathcal{A}$, and a sequence $x_{n}\subseteq\mathbb{R}^{n}$
such that $\Omega_{n}+x_{n}\to\Omega$ in the metric $d$.

Note that, since our energy is translation invariant, the above convergence
result is sufficient for our purposes.

###### Proof.

In this proof it is more convenient to work with constant speed, instead of
arc-length, parameterizations.

Consider minimizing sequence $\\{\Omega_{n}\\}\subseteq\mathcal{A}$, and let
$\varphi_{n}:[0,1]\longrightarrow\partial\Omega_{n}$ be constant speed
parameterizations. Note all $\partial\Omega_{n}$ are closed curves, and as
$E_{p,\lambda}$ is translation invariant, we can replace $\Omega_{n}$ with
translated copies (which, for brevity, we still denote by $\Omega_{n}$, and by
$\varphi_{n}$ the parameterization of $\partial\Omega_{n}$) such that
$\varphi_{n}(0)=\varphi_{n}(1)=0$. We show that we are under the conditions
(3.12): first, the upper bound on the perimeter (2.7) and
$\Omega_{n}\subseteq\mathbb{R}^{2}$ ensures all $\Omega_{n}$ are contained in
some compact set $K$. As the curves $\varphi_{n}$ are parameterized by
constant speed, we have
$\|\varphi_{n}^{\prime}\|=L(\varphi_{n})=\mathcal{H}^{1}(\partial\Omega_{n})$
a.e. Then, in view of Lemma 2.2, we infer (3.12). Therefore there exists a
limit curve $\varphi:[0,1]\longrightarrow K$ such that the convergences in
Lemma 3.1 hold. Since $\varphi_{n}(0)=\varphi_{n}(1)=0$ for all $n$, we get
$\varphi(0)=\varphi(1)=0$ too. We define $\Omega$ to be the bounded area
delimited by $\varphi$, and the graph of $\varphi$ turn out to be
$\partial\Omega$. By construction, $\Omega$ is compact.

We need to check it is convex: consider arbitrary $P,Q\in\Omega$, $t\in(0,1)$,
and we show that $(1-t)P+tQ\in\Omega$. Consider sequences
$P_{n},Q_{n}\in\Omega_{n}$ such that $P_{n}\to P$, $Q_{n}\to Q$: since each
$\Omega_{n}$ is convex, $(1-t)P_{n}+tQ_{n}\in\Omega_{n}$. By Lemma 3.1, we
know $\|\varphi_{n}-\varphi\|_{C^{0}([0,1];\mathbb{R}^{2})}\to 0$. As a
consequence,

$d_{\mathcal{H}}(\partial\Omega_{n},\partial\Omega)\to 0$

too, This allows us to choose, for each $n$, another point $z_{n}\in\Omega$
such that $|z_{n}-((1-t)P_{n}+tQ_{n})|\leq
d_{\mathcal{H}}(\partial\Omega_{n},\partial\Omega)$. By construction, now the
sequences $(1-t)P_{n}+tQ_{n}$ and $z_{n}$ have the same limit. As
$(1-t)P_{n}+tQ_{n}\to(1-t)P+tQ$, and $z_{n}\to z$, hence $z=(1-t)P+tQ$, using
the compactness of $\Omega$ finally gives $z\in\Omega$.

Finally, we check that $\dim_{\mathcal{H}}\Omega=2$. Since the ambient space
$\mathbb{R}^{2}$ has already Hausdorff dimension two, it suffices to show that
$\Omega$ contains a set of Hausdorff dimension two. For each $n$, we can use
the construction from the proof of Lemma 2.2 on each $\Omega_{n}$: we showed
the existence of triangles $T_{i}:=\triangle\tilde{A}\tilde{P}_{i}\tilde{B}$
(see Figure 1) whose distance to the boundary is at least $r_{i}/2$, with
$r_{i}$ being the incenter which satisfied $r_{i}\geq h_{i}/3$. Now, since we
showed in the proof of Lemma 2.2 that

$\sum_{i=1}^{2}h_{i,n}\operatorname{{\operatorname{diam}}}\Omega_{n}\geq\mathcal{H}^{2}(\Omega_{n}),$

and

$\mathcal{H}^{2}(\Omega_{n})\geq
C_{1},\qquad\operatorname{{\operatorname{diam}}}\Omega_{n}\leq\mathcal{H}^{1}(\partial\Omega_{n})\leq
C_{3}$

due to Corrollary 2.1, we get

$\sum_{i=1}^{2}h_{i,n}\geq\frac{\mathcal{H}^{2}(\Omega_{n})}{\operatorname{{\operatorname{diam}}}\Omega_{n}}\geq\frac{C_{1}}{C_{3}}>0.$

This shows that at least one of the triangles $T_{i,n}$, $i=1,2$, must be non
degenerate since its inradius is bounded from below by

$\max_{i=1,2}r_{i,n}\geq\max_{i=1,2}\frac{h_{i,n}}{3}\geq\frac{C_{1}}{6C_{3}},$

and the proof is complete. ∎

Lemma 3.3 is of crucial importance: since we are interested in the minimizers
of $E_{p,\lambda}$, this allows us to reduce the minimization problem to
$\mathcal{A}$, and neglect the highly irregular elements of
$\bar{\mathcal{A}}\setminus\mathcal{A}$.

######  Lemma 3.4.

Given $p\geq 1$, $\lambda>0$, and a minimizing sequence
$\Omega_{n}\subseteq{\mathcal{A}}$ converging to $\Omega\in{\mathcal{A}}$ with
respect to $d$, then it holds:

$\displaystyle\mathcal{H}^{2}(\Omega)$
$\displaystyle=\lim_{n\to+\infty}\mathcal{H}^{2}(\Omega_{n}),$ (3.13)
$\displaystyle\mathcal{H}^{1}(\partial\Omega)$
$\displaystyle\leq\liminf_{n\to+\infty}\mathcal{H}^{1}(\partial\Omega_{n}),$
(3.14)
$\displaystyle\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
$\displaystyle=\lim_{n\to+\infty}\int_{\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x.$
(3.15)

###### Proof.

Estimate (3.13) follows from the definition of the metric $d$ and Remark 3.2.

To prove (3.14), recall that the perimeter
$\mathcal{H}^{1}(\partial\Omega_{n})$ is the total variation of the
characteristic function of $\Omega_{n}$. Convergence $\Omega_{n}\to\Omega$
with respect to $d$ implies (see Remark 3.2)

$\chi_{\Omega_{n}}\to\chi_{\Omega}\text{ strongly in }L^{1}(\mathbb{R}^{2}),$

with “$\chi$” denoting the characteristic function of the subscribed set. Thus
(3.14) follows from the lower-semicontinuity of the total variation semi-norm.

To prove (3.15), note that

$\displaystyle\int_{\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x$
$\displaystyle=\int_{\Omega_{n}\backslash\Omega}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x+\int_{\Omega_{n}\cap\Omega}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x$
$\displaystyle\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
$\displaystyle=\int_{\Omega\backslash\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\int_{\Omega_{n}\cap\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x,$

hence

$\displaystyle\bigg{|}\int_{\Omega_{n}}$
$\displaystyle\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x-\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x\bigg{|}$
$\displaystyle\leq\int_{\Omega_{n}\backslash\Omega}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x+\int_{\Omega\backslash\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
(3.16)
$\displaystyle+\int_{\Omega_{n}\cap\Omega}|\operatorname{dist}^{p}(x,\partial\Omega_{n})-\operatorname{dist}^{p}(x,\partial\Omega)|\,{\operatorname{d}}x.$
(3.17)

By Lemma 2.2,

$\operatorname{{\operatorname{diam}}}(\Omega_{n})\leq\mathcal{H}^{1}(\partial\Omega_{n})\leq
C_{3}.$

According to (3.14),

$\operatorname{{\operatorname{diam}}}(\Omega)\leq\mathcal{H}^{1}(\partial\Omega)\leq\liminf_{n\to+\infty}\mathcal{H}^{1}(\partial\Omega_{n})\leq
C_{3}.$

Therefore,

$\displaystyle\int_{\Omega_{n}\backslash\Omega}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x$
$\displaystyle\leq\mathcal{H}^{2}(\Omega_{n}\backslash\Omega)(\operatorname{{\operatorname{diam}}}(\Omega_{n}))^{p}\leq\mathcal{H}^{2}(\Omega_{n}\backslash\Omega)C_{3}^{p}\to
0,$
$\displaystyle\int_{\Omega\backslash\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x$
$\displaystyle\leq\mathcal{H}^{2}(\Omega\backslash\Omega_{n})(\operatorname{{\operatorname{diam}}}(\Omega))^{p}\leq\mathcal{H}^{2}(\Omega\backslash\Omega_{n})C_{3}^{p}\to
0,$

hence the sum in (3.16) goes to zero. To estimate (3.17), denote by
$d_{\mathcal{H}}$ the Hausdorff distance, and note that, by the Mean Value
theorem, it holds

$\displaystyle\int_{\Omega_{n}\cap\Omega}$
$\displaystyle|\operatorname{dist}^{p}(x,\partial\Omega_{n})-\operatorname{dist}^{p}(x,\partial\Omega)|\,{\operatorname{d}}x$
$\displaystyle\leq\int_{\Omega_{n}\cap\Omega}|\operatorname{dist}(x,\partial\Omega_{n})-\operatorname{dist}(x,\partial\Omega)|$
$\displaystyle\cdot
p\sup_{x\in\Omega_{n}\cap\Omega}\Big{(}\max\\{\operatorname{dist}(x,\partial\Omega_{n}),\operatorname{dist}(x,\partial\Omega)\\}\Big{)}^{p-1}\,{\operatorname{d}}x$
$\displaystyle\leq\mathcal{H}^{2}(\Omega_{n}\cap\Omega)d_{\mathcal{H}}(\partial\Omega_{n},\partial\Omega)\cdot
p\Big{(}\max\\{\operatorname{{\operatorname{diam}}}\Omega_{n},\operatorname{{\operatorname{diam}}}\Omega\\}\Big{)}^{p-1}$
$\displaystyle\leq\mathcal{H}^{2}(\Omega_{n}\cap\Omega)\
d_{\mathcal{H}}(\partial\Omega_{n},\partial\Omega)\cdot p\ C_{3}^{p-1}\to 0.$

Thus the term in (3.17) goes to zero too, and (3.15) is proven. ∎

Now we prove part (1) of Theorem 1.1, i.e., the existence of minimizers in
${\mathcal{A}}$.

######  Lemma 3.5.

For any $p\geq 1$, $\lambda>0$, the functional $E_{p,\lambda}$ admits a
minimizer $\Omega\in{\mathcal{A}}$, which satisfies:

$\displaystyle\mathcal{H}^{2}(\Omega)\geq
C_{1},\qquad\mathcal{H}^{1}(\partial\Omega)$ $\displaystyle\leq C_{3},$

with $C_{1}$ (resp. $C_{3}$) defined in (2.3) (resp. (2.7)).

###### Proof.

Corollary 2.1 gives $\mathcal{H}^{2}(\Omega_{n})\geq C_{1}$ for any
sufficiently large $n$, and Lemma 3.4 gives

$\mathcal{H}^{2}(\Omega)=\lim_{n\to+\infty}\mathcal{H}^{2}(\Omega_{n})\geq
C_{1}.$ (3.18)

Based on Lemma 3.4 and equation (2.7),

$\displaystyle\mathcal{H}^{1}(\partial\Omega)\leq C_{3}.$ (3.19)

Lemma 3.4 gives

$\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\leq\liminf_{n\to+\infty}\frac{\mathcal{H}^{1}(\partial\Omega_{n})}{\mathcal{H}^{2}(\Omega_{n})}.$
(3.20)

and

$\displaystyle\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x=\lim_{n\to+\infty}\int_{\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x.$
(3.21)

Combining (3.20) and (3.21) gives

$\displaystyle E_{p,\lambda}(\Omega)$
$\displaystyle=\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}$
$\displaystyle\leq\lim_{n\to+\infty}\int_{\Omega_{n}}\operatorname{dist}^{p}(x,\partial\Omega_{n})\,{\operatorname{d}}x+\lambda\liminf_{n\to+\infty}\
\frac{\mathcal{H}^{1}(\partial\Omega_{n})}{\mathcal{H}^{2}(\Omega_{n})}\leq\liminf_{n\to+\infty}E_{p,\lambda}(\Omega_{n})=\inf_{\bar{\mathcal{A}}}E_{p,\lambda},$

hence $\Omega$ is effectively a minimizer of $E_{p,\lambda}$ in
$\bar{\mathcal{A}}$. Lemma 3.3 shows $\Omega\in{\mathcal{A}}$.

∎

##  4 Regularity

Now we prove part (2) of Theorem 1.1. The proof will be split over Lemmas 4.2
and 4.3.

######  Lemma 4.1.

Let $S$ be a compact, convex set, with Hausdorff dimension equal to 2. Let
$w_{1}$, $w_{2}\in\partial S$ be arbitrary distinct points, and let $\sigma$
be the segment with endpoints $w_{1}$ and $w_{2}$. Denoting by $S_{1}$ and
$S_{2}$ the two connected components of $S\backslash\sigma$, then both
$S_{1}$, $S_{2}$ are convex.

###### Proof.

Endow $\mathbb{R}^{2}$ with a Cartesian coordinate system. Upon rotation and
reflection, assume that $\sigma$ lies in the $y$-axis, and
$S_{1}\subseteq\\{x>0\\}$, $S_{2}\subseteq\\{x<0\\}$. Clearly, given points
$u,v\in S_{1}$, the segment $\xi$ between $u$ and $v$ lies entirely in
$S\cap\\{x>0\\}=S_{1}$, hence $S_{1}$ is convex. The proof for $S_{2}$ is
analogous. ∎

######  Lemma 4.2.

($C^{1}$-regularity) For any $p\geq 1$, $\lambda>0$, any minimizer of
$E_{p,\lambda}$ is $C^{1}$-regular.

###### Proof.

Consider an arbitrary minimizer $\Omega\in{\mathcal{A}}$. Endow
$\mathbb{R}^{2}$ with a polar coordinate system. We parameterize
$\partial\Omega$ by a closed Lipschitz curve

$\gamma:[0,2\pi]\longrightarrow\partial\Omega.$

The proof is achieved by a contradiction argument. Assume that $\Omega$ is not
$C^{1}$-regular. That is, $\gamma$ is not $C^{1}$-regular at some point
$t_{0}$. Upon rotating the coordinates, we can also assume $t_{0}\in(0,2\pi)$.
Since $\Omega$ is convex, both one-sided derivatives

$l^{-}:=\lim_{t\to t_{0}^{-}}\gamma^{\prime}(t),\qquad l^{+}:=\lim_{t\to
t_{0}^{+}}\gamma^{\prime}(t)$

are well-defined [1, 20]. Denote by $\alpha$ the angle between $l^{-}$ and
$l^{+}$. Clearly, $\alpha\neq\pi$.

Figure 3: A schematic representation (near $\gamma(t_{0})$, in first order
approximation in $\varepsilon$) of the construction of $\Omega_{\varepsilon}$.

Figure 3 is a representation (in first order approximation) of
$\partial\Omega$ near $\gamma(t_{0})$. For small parameters $0<\varepsilon\ll
1$, construct the competitor $\Omega_{\varepsilon}$ as follows:

  1. 1.

Choose $t_{1}<t_{0}<t_{2}$ such that (in first order approximation in
$\varepsilon$)

$\mathcal{H}^{1}(\gamma([t_{1},t_{0}]))=\mathcal{H}^{1}(\gamma([t_{0},t_{2}]))=\varepsilon+O(\varepsilon^{2}).$

  2. 2.

Denote by

$\sigma:=\\{(1-s)\gamma(t_{1})+s\gamma(t_{2}):s\in[0,1]\\}$

the line segment between $\gamma(t_{1})$ and $\gamma(t_{2})$, and set

$L:=\big{(}\partial\Omega\backslash\gamma([t_{1},t_{2}])\big{)}\cup\sigma.$
(4.22)

Note that such $L$ is a convex Jordan curve, and denote by
$\Omega_{\varepsilon}$ the bounded region delimited by $L$.

By construction, in first order approximation in $\varepsilon$, it holds

$\displaystyle\mathcal{H}^{1}(\partial\Omega_{\varepsilon})$
$\displaystyle=\mathcal{H}^{1}(\partial\Omega)-2\varepsilon(1-\sin(\alpha/2))+O(\varepsilon^{2}),$
(4.23) $\displaystyle\mathcal{H}^{2}(\Omega_{\varepsilon})$
$\displaystyle=\mathcal{H}^{2}(\Omega)-\frac{\varepsilon^{2}\sin\alpha}{2}+o(\varepsilon^{2})=\mathcal{H}^{2}(\Omega)+O(\varepsilon^{2}).$
(4.24)

Moreover, it is straightforward to show that

$\int_{\Omega_{\varepsilon}}\operatorname{dist}^{p}(x,\partial\Omega_{\varepsilon})\,{\operatorname{d}}x\leq\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x.$
(4.25)

Recalling that $\mathcal{H}^{2}(\Omega)>0$ (since $\Omega$ is a minimizer),
combining (4.23), (4.24) and (4.25) gives (in first order approximation in
$\varepsilon$)

$\displaystyle E_{p,\lambda}(\Omega_{\varepsilon})$
$\displaystyle=\int_{\Omega_{\varepsilon}}\operatorname{dist}^{p}(x,\partial\Omega_{\varepsilon})\
dx+\lambda\frac{\mathcal{H}^{1}(\partial\Omega_{\varepsilon})}{\mathcal{H}^{2}(\Omega_{\varepsilon})}$
$\displaystyle\leq\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\
dx+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)-2\varepsilon(1-\sin(\alpha/2))+O(\varepsilon^{2})}{\mathcal{H}^{2}(\Omega)+O(\varepsilon^{2})}$
$\displaystyle=\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}-\frac{2\lambda\varepsilon(1-\sin(\alpha/2))}{\mathcal{H}^{2}(\Omega)}+O(\varepsilon^{2})$
$\displaystyle=E_{p,\lambda}(\Omega)-\frac{2\lambda\varepsilon(1-\sin(\alpha/2))}{\mathcal{H}^{2}(\Omega)}+O(\varepsilon^{2})$
$\displaystyle=\min_{\mathcal{A}}E_{p,\lambda}-\frac{2\lambda\varepsilon(1-\sin(\alpha/2))}{\mathcal{H}^{2}(\Omega)}+O(\varepsilon^{2}),$

which is a contradiction for sufficiently small $\varepsilon$. Thus $\Omega$
must be $C^{1}$-regular. ∎

######  Lemma 4.3.

Given $p\geq 1$, $\lambda>0$, a minimizer $\Omega$ of $E_{p,\lambda}$, let
$\gamma:[0,\mathcal{H}^{1}(\partial\Omega)]\longrightarrow\partial\Omega$ be
an arc-length parameterization. Then it holds

$\displaystyle\limsup_{h\to 0}$
$\displaystyle\frac{|\gamma(t+2h)-2\gamma(t)+\gamma(t-2h)|}{h^{2}}\leq 4C,$
(4.26)

for any $t$, where $C$ is some constant depending only on $\lambda$ and $p$
(and independent of $\Omega$).

We remark that (4.26) implies $C^{1,1}$-regularity of $\partial\Omega$.

###### Proof.

Consider an arbitrary point $p_{0}\in\partial\Omega$. Since we proved that
$\Omega$ is $C^{1}$-regular, consider a (local) orthogonal coordinate system
with origin in $p_{0}$, and $x$-axis oriented along the tangent derivative (at
$p_{0}$), such that $\Omega$ is entirely contained in the half-plane $\\{y\geq
0\\}$. The boundary $\partial\Omega$ is thus (locally) the graph of some
nonnegative function $f$. Clearly, such $f$ satisfies $f(0)=0$.

Figure 4: A schematic representation of the construction near $p_{0}=(0,0)$.

Choose an arbitrary $0<\varepsilon\ll 1$. Denote by

$\sigma_{\varepsilon}:=\\{(x,y):0\leq x\leq\varepsilon,\ y=x\cdot
f(\varepsilon)/\varepsilon\\}$

the segment between the origin and $(\varepsilon,f(\varepsilon))$. Let
$L_{\varepsilon}$ be the curve obtained by replacing $f([0,\varepsilon])$ with
$\sigma_{\varepsilon}$. That is,

$L_{\varepsilon}:=(\partial\Omega\backslash
f([0,\varepsilon]))\cup\sigma_{\varepsilon}.$

By construction (see Lemma 4.1) $L_{\varepsilon}$ is a convex Jordan curve,
and let $\Omega_{\varepsilon}$ be the bounded region delimited by
$L_{\varepsilon}$. Note that:

  1. 1.

Clearly we can infer

$\int_{\Omega_{\varepsilon}}\operatorname{dist}^{p}(x,\partial\Omega_{\varepsilon})\,{\operatorname{d}}x\leq\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x.$
(4.27)

  2. 2.

For areas, since by construction it holds
$\Omega_{\varepsilon}\subseteq\Omega$, we have

$\displaystyle\mathcal{H}^{2}(\Omega)-\mathcal{H}^{2}(\Omega_{\varepsilon})$
$\displaystyle=\mathcal{H}^{2}(\Omega\backslash\Omega_{\varepsilon})$
$\displaystyle=\mathcal{H}^{2}(\\{(x,y):0\leq x\leq\varepsilon,\ f(x)\leq
y\leq x\cdot f(\varepsilon)/\varepsilon\\})$
$\displaystyle=\int_{0}^{\varepsilon}[xf(\varepsilon)/\varepsilon-f(x)]\,{\operatorname{d}}x=\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x.$
(4.28)

  3. 3.

For perimeters, note that $f^{\prime}(0)=0$, so $|f^{\prime}|$ is small near
$0$. In particular, by choosing sufficiently small $\varepsilon$, we can
ensure that

$\sqrt{1+|f^{\prime}(x)|^{2}}\leq 2\qquad\text{for all }x\in(0,\varepsilon),$

and also $|f(\varepsilon)|/\varepsilon$ can be made as small as we need, so to
satisfy

$\sqrt{1+\frac{f(\varepsilon)^{2}}{\varepsilon^{2}}}=1+\frac{f(\varepsilon)^{2}}{2\varepsilon^{2}}-\frac{1}{8}\Big{(}\frac{f(\varepsilon)^{2}}{\varepsilon^{2}}\Big{)}^{2}+O\Big{(}\Big{(}\frac{f(\varepsilon)^{2}}{\varepsilon^{2}}\Big{)}^{3}\Big{)}\leq
1+\frac{f(\varepsilon)^{2}}{2\varepsilon^{2}}.$

Therefore,

$\displaystyle\mathcal{H}^{1}(\partial\Omega)-\mathcal{H}^{1}(\partial\Omega_{\varepsilon})$
$\displaystyle=\int_{0}^{\varepsilon}\Big{(}\sqrt{1+|f^{\prime}(x)|^{2}}-\sqrt{1+\frac{f(\varepsilon)^{2}}{\varepsilon^{2}}}\Big{)}\,{\operatorname{d}}x$
$\displaystyle\geq\int_{0}^{\varepsilon}\Big{(}\sqrt{1+|f^{\prime}(x)|^{2}}-1-\frac{f(\varepsilon)^{2}}{2\varepsilon^{2}}\Big{)}\,{\operatorname{d}}x,$

where, since for sufficiently small $\varepsilon\ll 1$ the quantity
$\frac{f(\varepsilon)^{2}}{2\varepsilon^{2}}$ can be made arbitrarily small,
we have

$\int_{0}^{\varepsilon}\frac{f(\varepsilon)^{2}}{2\varepsilon^{2}}\,{\operatorname{d}}x=o(\varepsilon),\qquad\varepsilon\ll
1.$

Thus, for all sufficiently small $\varepsilon$,

$\displaystyle\mathcal{H}^{1}(\partial\Omega)-\mathcal{H}^{1}(\partial\Omega_{\varepsilon})$
$\displaystyle=\int_{0}^{\varepsilon}\Big{(}\sqrt{1+|f^{\prime}(x)|^{2}}-\sqrt{1+\frac{f(\varepsilon)^{2}}{\varepsilon^{2}}}\Big{)}\,{\operatorname{d}}x$
$\displaystyle\geq\int_{0}^{\varepsilon}\Big{(}\sqrt{1+|f^{\prime}(x)|^{2}}-1\Big{)}\,{\operatorname{d}}x+o(\varepsilon)$
$\displaystyle=\int_{0}^{\varepsilon}\frac{|f^{\prime}(x)|^{2}}{\sqrt{1+|f^{\prime}(x)|^{2}}+1}\,{\operatorname{d}}x+o(\varepsilon)\geq\frac{1}{3}\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x.$
(4.29)

Combining (4.27), (4.28) and (4.29) gives

$\displaystyle E_{p,\lambda}(\Omega_{\varepsilon})$
$\displaystyle=\int_{\Omega_{\varepsilon}}\operatorname{dist}^{p}(x,\partial\Omega_{\varepsilon})\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial\Omega_{\varepsilon})}{\mathcal{H}^{2}(\Omega_{\varepsilon})}$
$\displaystyle\leq\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)-\frac{1}{3}\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}.$
(4.30)

Since

$\displaystyle\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}$
$\displaystyle=\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\cdot\frac{\mathcal{H}^{2}(\Omega)}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}$
$\displaystyle=\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\cdot\bigg{(}1+\frac{\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}\bigg{)},$

estimate (4.30) reads

$\displaystyle E_{p,\lambda}(\Omega_{\varepsilon})$
$\displaystyle\leq\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x-\lambda\frac{\frac{1}{3}\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}$
$\displaystyle+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\cdot\bigg{(}1+\frac{\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}\bigg{)}$
$\displaystyle=E_{p,\lambda}(\Omega)+\lambda\frac{\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)-\frac{1}{3}\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x}{\mathcal{H}^{2}(\Omega)-(\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x)}.$
(4.31)

Since $\Omega$ is a minimizer, Lemma 3.5 gives $\mathcal{H}^{2}(\Omega)>0$,
and note that

$\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x\leq\frac{\mathcal{H}^{2}(\Omega)}{2}$

for all sufficiently small $\varepsilon$, hence the denominator in (4.31) is
positive. Thus the minimality of $\Omega$ forces the numerator in (4.31) to be
nonnegative, i.e.,

$3\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\bigg{(}\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x\bigg{)}-\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x\geq
0.$ (4.32)

Equation (4.28) shows

$\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x=\mathcal{H}^{2}(\Omega)-\mathcal{H}^{2}(\Omega_{\varepsilon})\geq
0$

since by construction $\Omega_{\varepsilon}\subseteq\Omega$. Lemma 3.5 gives

$\mathcal{H}^{2}(\Omega)\geq C_{1},\qquad\mathcal{H}^{1}(\partial\Omega)\leq
C_{3},$

hence

$3\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}\leq\frac{3C_{3}}{C_{1}}=:C,$

and (4.32) forces

$C\bigg{(}\frac{f(\varepsilon)\varepsilon}{2}-\int_{0}^{\varepsilon}f(x)\,{\operatorname{d}}x\bigg{)}\geq\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x.$
(4.33)

Since $\Omega$ is convex, and we assumed (at the beginning of this proof) that
$\Omega\subseteq\\{y\geq 0\\}$, $f$ is nonnegative, hence (4.33) forces

$\frac{C}{2}f(\varepsilon)\varepsilon\geq\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x.$
(4.34)

Note that, since $f(0)=0$, it follows

$f(\varepsilon)=\int_{0}^{\varepsilon}f^{\prime}(x)\,{\operatorname{d}}x\leq\int_{0}^{\varepsilon}|f^{\prime}(x)|\,{\operatorname{d}}x.$
(4.35)

By Hölder’s inequality,

$\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x\geq\frac{1}{\varepsilon}\bigg{(}\int_{0}^{\varepsilon}|f^{\prime}(x)|\,{\operatorname{d}}x\bigg{)}^{2},$
(4.36)

hence

$\displaystyle\frac{C}{2}$
$\displaystyle\varepsilon\int_{0}^{\varepsilon}|f^{\prime}(x)|\,{\operatorname{d}}x\overset{\eqref{f1}}{\geq}\frac{C}{2}f(\varepsilon)\varepsilon\overset{\eqref{E4}}{\geq}\int_{0}^{\varepsilon}|f^{\prime}(x)|^{2}\,{\operatorname{d}}x\overset{\eqref{fH}}{\geq}\frac{1}{\varepsilon}\bigg{(}\int_{0}^{\varepsilon}|f^{\prime}(x)|\,{\operatorname{d}}x\bigg{)}^{2}$
$\displaystyle\Longrightarrow\frac{C}{2}\varepsilon^{2}\geq\int_{0}^{\varepsilon}|f^{\prime}(x)|\,{\operatorname{d}}x\geq\int_{0}^{\varepsilon}f^{\prime}(x)\,{\operatorname{d}}x=f(\varepsilon).$

The above arguments can be repeated for $\varepsilon<0$, $|\varepsilon|\ll 1$
(or equivalently, when the orientation of $x$-axis is inverted). The
arbitrariness of $\varepsilon$ then gives

$\displaystyle\limsup_{\varepsilon\to 0}$
$\displaystyle\frac{|f(\varepsilon)-2f(0)+f(-\varepsilon)|}{(\varepsilon/2)^{2}}\leq
4C,$

concluding the proof. ∎

Now we prove part (3) of Theorem 1.1.

######  Lemma 4.4.

Given $p\geq 1$, $\lambda>0$, any minimizer $\Omega$ of $E_{p,\lambda}$
satisfies

$\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}=\frac{p+2}{\lambda(p+3)}\min_{\mathcal{A}}E_{p,\lambda}.$

###### Proof.

Let $\Omega$ be an arbitrary minimizer. Endow $\mathbb{R}^{2}$ with a
Cartesian coordinate system, and assume without loss of generality that
$(0,0)$ is in the interior part of $\Omega$. For any $r>0$, denote by

$T_{r}:\mathbb{R}^{2}\longrightarrow\mathbb{R}^{2},\qquad T_{r}(x):=rx$

the homothety of center $(0,0)$ and ratio $r$. Note that
$T_{r}(\Omega)\in{\mathcal{A}}$ for any $r>0$, and the scalings are

$\displaystyle\int_{T_{r}(\Omega)}\operatorname{dist}^{p}(x,\partial
T_{r}(\Omega))\,{\operatorname{d}}x$
$\displaystyle=r^{p+2}\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x,$
$\displaystyle\frac{\mathcal{H}^{1}(\partial
T_{r}(\Omega))}{\mathcal{H}^{2}(T_{r}(\Omega))}$
$\displaystyle=\frac{1}{r}\cdot\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}.$

Define the function

$f:(0,+\infty)\longrightarrow(0,+\infty),\qquad
f(r):=E_{p,\lambda}(T_{r}(\Omega))=r^{p+2}\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\frac{\lambda}{r}\cdot\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}.$

Since $f$ is smooth, and attains a global minimum at $r=1$, it follows

$\displaystyle f^{\prime}(1)$
$\displaystyle=(p+2)\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x-\lambda\cdot\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}=0$
$\displaystyle\Longrightarrow\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x=\frac{\lambda}{p+2}\cdot\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)},$

hence

$E_{p,\lambda}(\Omega)=\frac{\lambda(p+3)}{p+2}\cdot\frac{\mathcal{H}^{1}(\partial\Omega)}{\mathcal{H}^{2}(\Omega)}=\min_{\mathcal{A}}E_{p,\lambda},$

and the proof is complete. ∎

Let us conclude the paper with some final remarks. In this paper we
investigated the minimization problem for the average distance functional,
with perimeter-to-area ratio penalization, in the plane. We proved the
existence and $C^{1,1}$ regularity of minimizers, mainly relying on
constructing suitable competitors. Echoing and developing former studies that
exclusively focused on either the 1D average distance problem or purely
surface area-to-volume ratio question, by considering optimal sets of combined
energy from broader and more eclectic perspectives, this study enriches and
deepens our understanding of penalized average distance problem.

We remark that all the main results of this paper, i.e. bounds on the
perimeter and area, and $C^{1,1}$-regular of minimizers, can be also proven if
we replace the perimeter-to-area term with a generalized ratio of the form
$\lambda\frac{\mathcal{H}^{1}(\partial\Omega)^{\alpha}}{\mathcal{H}^{2}(\Omega)^{\beta}}$,
symbolizing a perimeter term normalized (by area) with different scaling
exponents $\alpha$ and $\beta$. That is, we consider an energy of the form

$E_{p,\lambda}^{\alpha,\beta}(\Omega):=\int_{\Omega}\operatorname{dist}^{p}(x,\partial\Omega)\,{\operatorname{d}}x+\lambda\frac{\mathcal{H}^{1}(\partial\Omega)^{\alpha}}{\mathcal{H}^{2}(\Omega)^{\beta}},$
(4.37)

where $\alpha,\beta$ are given powers satisfying
$2\beta>\alpha>\frac{p}{p+1}\beta>0$. This last bound, combined with Young’s
inequality, allows us to easily bound the perimeter, and the subsequent
results. It can also be quickly checked that if $\alpha>2\beta$, then
minimizers are just single points. One more remark is that according to (2.8),
if in (4.37) we pick $\alpha=p,\beta=p+1$ and $\lambda=C$ as in (2.8), we get

$\displaystyle E_{p,\lambda}^{p,p+1}(\Omega)\geq
C\frac{\mathcal{H}^{2}(\Omega)^{p+1}}{\mathcal{H}^{1}(\partial\Omega)^{p}}+C\frac{\mathcal{H}^{1}(\partial\Omega)^{p}}{\mathcal{H}^{2}(\Omega)^{p+1}}\geq
2C.$

So in this case if the optimal constant in (2.8) is obtained by a circle, the
optimal shape for (1.1) is a circle. An interesting question worthy further
consideration is if the circle would be the minimizer for other parameters, as
in similar discussions given in [21, 16, 15, 13]. Another natural question is
to ask if in general one may improve the $C^{1,1}$ regularity by combining the
established results with elliptic regularity theory, given that the variation
of the perimeter-to-area ratio leads to a system of second order differential
equations of the boundary parametrization.

In addition, it is interesting to improve the results of this paper to higher
dimensions, again with a generalized ratio penalization. However, the
geometric complexity of higher dimensional objects can increase significantly,
and more work is required to exclude more complicated sets (e.g.,
“tentacles”), which were not an issue in the planar case, thus we expected to
rely on rather different tools and arguments.

## Acknowledgments

The work of QD is supported in part by the National Science Foundation under
award DMS-2012562 and the ARO MURI Grant W911NF-15-1-0562. XYL acknowledges
the support of his NSERC Discovery Grant “Regularity of minimizers and pattern
formation in geometric minimization problems”.

## References

  * [1] L. Ambrosio, N. Fusco and D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, The Clarendon Press, Oxford University Press, New York, 2000.
  * [2] G. Buttazzo, E. Mainini, E. Stepanov, Stationary configurations for the average distance functional and related problems, Control Cybernet, 38 (2009), pp. 1107-1130.
  * [3] G. Buttazzo, E. Oudet, E. Stepanov, Optimal transportation problems with free Dirichlet regions, in: Variational methods for discontinuous structures, in: Progr. Nonlinear Differential Equations Appl., vol. 51, Birkhäuser, Basel, 2002, pp. 41-65.
  * [4] G. Buttazzo, A. Pratelli, S. Solimini, E. Stepanov, Optimal urban networks via mass transportation, Lecture Notes in Mathematics, Springer-Verlag, Berlin, 2009.
  * [5] G. Buttazzo, A. Pratelli, E. Stepanov, Optimal pricing policies for public transportation networks, SIAM J. Optim., 16 (2006), pp. 826-853.
  * [6] G. Buttazzo, F. Santambrogio, A model for the optimal planning of an urban area, SIAM J. Math. Anal., 37 (2005), pp. 514-530.
  * [7] G. Buttazzo, F. Santambrogio, A mass transportation model for the optimal planning of an urban region, SIAM Rev., 51(2009), pp. 593-610.
  * [8] G. Buttazzo, E. Stepanov, Optimal transportation networks as free Dirichlet regions for the Monge-Kantorovich problem, Ann. Sc. Norm. Super. Pisa Cl. Sci. (5), 2(4) (2003), pp. 631-678.
  * [9] G. Buttazzo, E. Stepanov, Minimization problems for average distance functionals, in: Calculus of variations: topics from the mathematical heritage of E. De Giorgi, vol. 14 of Quad. Mat., Dept. Math., Seconda Univ. Napoli, Caserta, 2004, pp. 48-83.
  * [10] N. Campbell, et al., Biology: Concepts & Connections, Benjamin Cummings, San Francisco, CA, 2003.
  * [11] M. P. do Carmo, Differential Geometry of Curves and Surfaces, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1976\.
  * [12] R. Eckhoff, Dust Explosions in the Process Industries (Third Edition), Gulf Professional Publishing, 2003.
  * [13] A. Figalli, N. Fusco, F. Maggi, V. Millot, M. Morini, Isoperimetry and Stability Properties of Balls with Respect to Nonlocal Energies, Commun. Math. Phys., 336(2015), 441-507.
  * [14] W. Gao, et al., Flame propagation mechanisms in dust explosions, J. Loss Prev. Process Ind., 36 (2015), pp. 186-194.
  * [15] M. Goldman, M. Novaga, M. Röger, Quantitative estimates for bending energies and applications to non-local variational problems, Proc. R. Soc. Edinb. A, 150-1 (2020), pp. 131-169.
  * [16] M. Goldman, M. Novaga, B. Ruffini, On minimizers of an isoperimetric problem with long-range interactions under a convexity constraint, Anal. PDE, 11-5 (2018), pp. 113-1142.
  * [17] W. Helfrich, Elastic properties of lipid bilayers: theory and possible experiments, Zeitschrift für Naturforschung. Teil C: Biochemie, Biophysik, Biologie, Virologie, 28 (1973), pp. 693-703.
  * [18] A. Lemenant, A presentation of the average distance minimizing problem, Zap. Nauchn. Sem. S.-Peterburg. Otdel. Mat. Inst. Steklov. (POMI), 390 (2011), pp. 117-146, 308.
  * [19] A. Lemenant, E. Mainini, On convex sets that minimize the average distance, ESAIM Control Optim. Calc. Var., 18 (2012), pp. 1049-1072.
  * [20] X. Lu, D. Slepčev, Properties of Minimizers of Average-Distance Problem via Discrete Approximation of Measures, SIAM J. Math. Anal., 45-5 (2013), pp. 3114-3131.
  * [21] C. Muratov, H. Knüpfer, On an Isoperimetric Problem with a Competing Nonlocal Term II: The General Case, Commun. Pure Appl. Math., 67-12 (2014), pp. 1974-1994.
  * [22] E. Paolini, E. Stepanov, Qualitative properties of maximum distance minimizers and average distance minimizers in ${\mathbb{R}}^{n}$, J. Math. Sci. (N.Y.), 122 (3) (2004), pp. 3290-3309.
  * [23] F. Santambrogio, P. Tilli,Blow-up of optimal sets in the irrigation problem, J. Geom. Anal., 15 (2005), pp. 343-362.
  * [24] K. Schmidt-Nielsen, Scaling: Why is Animal Size so Important? Cambridge University Press, New York, NY, 1984.
  * [25] D. Slepčev, Counterexample to regularity in average-distance problem, Ann. Inst. H. Poincaré Anal. Non Linéaire, 31 (2014), pp. 169-184.
  * [26] P. Tilli, Some explicit examples of minimizers for the irrigation problem, J. Convex Anal., 17 (2010), pp. 583-595.
  * [27] S. Vogel, Life’s Devices: The Physical World of Animals and Plants, Princeton University Press, Princeton, NJ, 1988.
  * [28] R.A. Wijsman, Convergence of Sequences of Convex Sets, Cones and Functions. II, Trans. Amer. Math. Soc., 123-1 (1966), pp. 32-45.