UNESCO-UNISA-iTLABS Africa Chair in Nanoscience and Nanotechnology (U2ACN2),
College of Graduate Studies, University of South Africa (UNISA), Pretoria,
South Africa

# Screening of BindingDB database ligands against EGFR, HER2, Estrogen,
Progesterone and NF-$\kappa$B receptors based on machine learning and
molecular docking

Parham Rezaee  Department of Biophysics, School of Biological Sciences,
Tarbiat Modares University, Tehran, Iran  Shahab Rezaee  Department of
Biophysics, School of Biological Sciences, Tarbiat Modares University, Tehran,
Iran  Malik Maaza  UNESCO-UNISA-iTLABS Africa Chair in Nanoscience and
Nanotechnology (U2ACN2), College of Graduate Studies, University of South
Africa (UNISA), Pretoria, South Africa  Seyed Shahriar Arab  Department of
Pediatrics, University of California, San Diego, La Jolla, CA 92093, USA
<EMAIL_ADDRESS>
###### Abstract

Breast cancer, the second most prevalent cancer among women worldwide,
necessitates the exploration of novel therapeutic approaches. To target the
four subgroups of breast cancer “hormone receptor-positive and HER2-negative,
hormone receptor-positive and HER2-positive, hormone receptor-negative and
HER2-positive, and hormone receptor-negative and HER2-negative” it is crucial
to inhibit specific targets such as EGFR, HER2, ER, NF-$\kappa$B, and PR.

In this study, we evaluated various methods for binary and multiclass
classification. Among them, the GA-SVM-SVM:GA-SVM-SVM model was selected with
an accuracy of 0.74, an F1-score of 0.73, and an AUC of 0.94 for virtual
screening of ligands from the BindingDB database. This model successfully
identified 4454, 803, 438, and 378 ligands with over 90% precision in both
active/inactive and target prediction for the classes of EGFR+HER2, ER,
NF-$\kappa$B, and PR, respectively, from the BindingDB database. Based on to
the selected ligands, we created a dendrogram that categorizes different
ligands based on their targets. This dendrogram aims to facilitate the
exploration of chemical space for various therapeutic targets.

Ligands that surpassed a 90% threshold in the product of activity probability
and correct target selection probability were chosen for further investigation
using molecular docking. The binding energy range for these ligands against
their respective targets was calculated to be between -15 and -5
$\mathrm{kcal}\text{/}\mathrm{mol}$. Finally, based on general and common
rules in medicinal chemistry, we selected 2, 3, 3, and 8 new ligands with high
priority for further studies in the EGFR+HER2, ER, NF-$\kappa$B, and PR
classes, respectively.

###### keywords:

virtual screening, machine learning, molecular docking, breast cancer

##  1 Introduction

Breast cancer, characterized by the highest mortality rate among various
cancer types, is a widespread condition. The development and progression of
breast cancer are facilitated by the interaction of estrogen and progesterone
receptors with breast cells1. These hormones, estrogen and progesterone, bind
to their respective receptors in the cytoplasm, leading to dimerization and
subsequent entry into the nucleus. Additionally, they bind to estrogen and
progesterone response elements located near the promoters of target genes. In
a study by Shirazi, it was demonstrated that estradiol, an estrogen hormone,
alone promoted the growth of MCF-7 cells compared to the control group.
Conversely, tamoxifen, a well-known estrogen blocker, arrested the
proliferation of MCF-7 cells for a minimum of five days. When both the
stimulator (estradiol) and blocker (tamoxifen) were applied together, the
level of stimulation in MCF-7 cell growth was reduced. Consequently, both
estrogen and progesterone contribute to the initiation and acceleration of
breast cancer. An excessive presence of these receptors indicates hormone
(ER+/PR+) mediated breast cancer. Targeting estrogen and progesterone hormones
can aid in the identification of potent inhibitors for breast cancer2, 3.

Figure 1: The flowchart of these procedures.

EGFR is a transmembrane glycoprotein with a crucial role in cell signaling. It
consists of a ligand-binding domain and a tyrosine kinase domain. Binding with
endogenous epidermal growth factor (EGF) leads to dimerization,
autophosphorylation, and activation of downstream pathways, resulting in cell
proliferation and differentiation4. EGFR is a major target in breast cancer,
and anti-EGFR agents have shown efficacy, especially in patients with specific
EGFR mutations. Tyrosine kinase inhibitors (TKIs) like gefitinib, erlotinib,
and lapatinib have been used to inhibit EGFR overexpression5, 6.

HER2 is a protein with tyrosine kinase activity encoded by the
neu/erbB2/c-erbB2 oncogene. It belongs to the EGFR family and plays a crucial
role in the development of normal and malignant breast tissue. Approximately
30% of human breast carcinomas exhibit HER2 amplification7. HER2 interacts
with insulin-like growth factor receptor-1 and estrogen receptor, initiating
cell signaling. Targeting HER2 has been a focus in anti-cancer drug
development, with the discovery of inhibitors like neratinib and afatinib.
However, further trials are needed to confirm their efficacy against breast
cancer. HER2 overexpression can increase the activity of MMP-2 and MMP-9
proteases, promoting invasiveness of breast cancer cells8. It also amplifies
VEGF expression and MMP-9 activity, potentially triggering angiogenic
responses. Clinical studies have shown that breast cancer cases with Erbb2
gene amplification have reduced responsiveness to certain treatments compared
to cases with normal ErbB2 expression9, 10.

Nuclear factor-kappa B (NF-$\kappa$B) is a transcription factor that regulates
the expression of genes involved in cell proliferation, immunological
responses, and inflammation. It contributes to the development of breast
tumors, lymphoma, and colorectal cancer11. In breast cancer, NF-$\kappa$B
activation occurs downstream of EGFR signaling, specifically in the ER-
negative subtype. HER-2 overexpression leads to the activation of the PI3K/Akt
pathway and induction of NF-$\kappa$B. NF-$\kappa$B plays a role in
angiogenesis by stimulating the expression of VEGF and IL-8. It can activate
two signaling pathways: the classical (canonical) pathway and the alternative
(noncanonical) pathway12. Drugs like lapatinib and certain microtubule
disruptors activate NF-$\kappa$B. In vitro studies have shown that ginseng
inhibits COX-2 and NF-$\kappa$B activation in breast cancer cell lines13.

Despite remarkable advancements in the field of basic life sciences and
biotechnology, the process of drug discovery and development (DDD) continues
to be slow and costly. On average, it takes around 15 years and approximately
US$2 billion to develop a small-molecule drug14. While clinical studies are
widely acknowledged as the most expensive phase in drug development, the
greatest potential for time and cost savings lies in the earlier stages of
discovery and preclinical research. Preclinical efforts alone account for over
43% of pharmaceutical expenses, in addition to significant public funding14,
15. This is primarily due to the high attrition rate observed at every step,
ranging from target selection and hit identification to lead optimization and
the selection of clinical candidates. Furthermore, the substantial failure
rate in clinical trials, currently at 90%16, can largely be attributed to
issues originating in the early stages of discovery, such as inadequate target
validation or suboptimal properties of ligands. Discovering faster and more
accessible methods to identify a broader range of high-quality chemical
probes, hits, and leads with optimal absorption, distribution, metabolism,
excretion, and toxicology (ADMET) as well as pharmacokinetics (PK) profiles
during the early phases of DDD would significantly enhance outcomes in
preclinical and clinical studies. Consequently, this would enable the
development of more effective, accessible, and safer drugs17, 18.

In this study, first, we developed two models with the purpose of identifying
active/inactive molecules and determining the target of each ligand.
Subsequently, we constructed a pipeline using these models. This pipeline was
utilized to screen the bindingDB ligands, employing various thresholds of
model precision. The selected ligands were subjected to molecular docking to
assess their binding energy with their respective targets. Additionally, we
applied several established principles in medicinal chemistry to prioritize
further investigation of these selected molecules, such as molecular dynamics,
in vitro and in vivo studies. Furthermore, we examined the significance of the
features employed in creating the target predictor model, aiming to identify a
simple rule for acceptable accurate target recognition as a common rule. The
procedure’s flowchart of this study is illustrated in Figure 1.

##  2 Materials and methods

In this study, our dataset consisted of inhibitors targeting various breast
cancer targets: 7341 for EGFR, 2182 for HER2, 1859 for ER, 1273 for
NF-$\kappa$B, and 1439 for PR. To obtain these compounds, we downloaded five
specific sdf files from the Binding database19 website, each file
corresponding to a particular class of inhibitors. These sdf files were then
converted to gjf files using the OpenBabel20 software. The 3D structures of
all the molecules were optimized using the Austin model 1 Hamiltonian
implemented in Gaussian software21. After optimization, the molecules were
used to calculate molecular descriptors with the help of Alvadesc22 software.
A total of 5668 descriptors, including 0-, 1-, 2-, and 3D descriptors, were
generated. To streamline the dataset, descriptors with constant values in 90%
of the compounds were removed. Additionally, among descriptors with a
correlation above 0.9, the one exhibiting higher pair correlation with all
other descriptors was kept and the others were automatically excluded.
Following these processes, 1461 descriptors remained for further analysis.

Table 1: Number of active and inactive molecules for each class Target | active | inactive | total
---|---|---|---
EGFR | 4922 | 2419 | 7341
ER | 1223 | 636 | 1859
HER2 | 1393 | 789 | 2182
NF-$\kappa$B | 447 | 826 | 1273
PR | 1055 | 384 | 1439
total | 9040 | 5054 | 14094

Each sdf file contains activity information pertaining to a specific molecule,
indicating the affinity of that molecule towards different therapeutic
targets. We extracted the activity information from the downloaded sdf files
for each class of molecules, and saved it in separate vectors. The enumeration
of the collected data can be found in Table 1. Molecules with $IC_{50}$,
$K_{i}$, and $EC_{50}$ values below 2000 nM were categorized as active
inhibitors, while those with values exceeding 2000 nM were considered
inactive. Also we removed molecules with activity more than 10000 nM as
outlier data from the inactive dataset. Both active and inactive molecules
were utilized to develop and evaluate the active/inactive classifiers. These
models serve the purpose of screening extensive databases and identifying new
potent molecules for the treatment of breast cancer.

A dataset of 1039519 molecules was gathered from the BindingDB database. The
same preparation process used for the breast cancer inhibitors described
earlier was applied to these downloaded molecules. Additionally, the
descriptors selected for the breast cancer inhibitors were also chosen for
these molecules. This resulted in a data matrix of size $1039519\times 1461$,
which was used for further analysis.

We employed various methods, including k-best, K-Nearest Neighbors (KNN),
Gaussian Naive Bayes (GNB), Quadratic Discriminant Analysis (QDA), Random
Forest (RF), and Support Vector Machine (SVM), to independently select
descriptors for the active/inactive and target classifiers. To optimize the
selection of descriptors across all data, we utilized a Genetic Algorithm (GA)
in an optimal manner. The GA started with a population size of 200 and evolved
through a maximum of 1000 generations, employing a crossover rate of 0.5 and a
mutation rate of 0.2. The estimator was configured with the aforementioned
methods, utilizing 5-fold cross-validation and an accuracy scoring function.
The only difference between the feature selection processes of the
active/inactive and target classifiers was the maximum number of features. The
active/inactive classifier allowed a maximum of 64 features, while the target
classifier allowed a maximum of 128 features.

We utilized the chosen features to create an optimized binary classifier for
predicting active/inactive molecules. Various methods, including K-nearest
neighbors (KNN), support vector machine (SVM), decision tree (DT), random
forest (RF), naive Bayes (NB), linear discriminant analysis (LDA), and
quadratic discriminant analysis (QDA), were employed for this purpose. To
maximize the performance of each method, we conducted a grid search to
identify the best parameters for constructing the model. Since the data
reality is not imbalanced and due to a lack of studies on active molecules, we
selected 461 active and 461 inactive molecules from the dataset to create a
balanced training dataset for binary classification.

For the target classifier, we employed methods such as KNN, SVM, DT, logistic
regression (LR), RF, NB, Gaussian naive Bayes (GNB), LDA, and QDA. Similar to
the binary classifier, we utilized grid search to identify the optimal
parameters for constructing the model. To ensure balance in the training
dataset, we selected 132 active molecules for each class for the multiple
classifier as well.

Table 2: The evaluation of GA-SVM-SVM and GA-RF-RF for binary and GA-SVM-SVM and GA-QDA-SVM models for therapeutic classification. GA-SVM-SVM | precision | recall | f1-score | support
---|---|---|---|---
active | 0.78 | 0.70 | 0.74 | 461
inactive | 0.73 | 0.80 | 0.76 | 461
GA-RF-RF | precision | recall | f1-score | support
active | 0.75 | 0.73 | 0.74 | 461
inactive | 0.74 | 0.76 | 0.75 | 461
GA-SVM-SVM | precision | recall | f1-score | support
EGFR+HER2 | 0.95 | 0.92 | 0.93 | 132
ER | 0.91 | 0.95 | 0.93 | 132
NF-$\kappa$B | 0.94 | 0.96 | 0.95 | 132
PR | 0.96 | 0.93 | 0.95 | 132
GA-QDA-SVM | precision | recall | f1-score | support
EGFR+HER2 | 0.94 | 0.92 | 0.93 | 132
ER | 0.95 | 0.95 | 0.95 | 132
NF-$\kappa$B | 0.93 | 0.94 | 0.93 | 132
PR | 0.96 | 0.96 | 0.96 | 132

Based on the selected best models for the active/inactive and target
classifiers, they were combined into a single pipeline as a decision model for
predicting the activity and target of each BindingDB ligand. The decision-
making process for the model’s predictions of activity and target was
constrained by a certainty threshold of 0.8, 0.85, and 0.9 for the
active/inactive classifier, and 0.9 for the target classifier.

Autodock Vina23 was utilized for molecular docking to calculate the binding
affinities between ligands and their respective targets. The grid box
resolution was set with specific coordinates for each target: (EGFR)(PDB ID:
1M17) had coordinates of 23.424, 1.310, 51.002 along the x, y, and z axes,
respectively, with a grid spacing of 0.2
$\mathrm{\SIUnitSymbolAngstrom}\text{/}$; (HER2)(PDB ID: 3PP0) had coordinates
of 17.563, 16.689, 26.321; (PR)(PDB ID: 1A28) had coordinates of 17.038,
0.145, 74.798; (ER)(PDB ID: 2IOK) had coordinates of 19.050, 35.696, 52.244;
and (NF-$\kappa$B)(PDB ID: 4KIK) had coordinates of 48.268, 31.589, -57.885.
These coordinates were used to define the binding sites for the docking
process. The grid dimensions were set at $25.2\times 25.2\times 25.2$
$\mathrm{\SIUnitSymbolAngstrom}\text{/}$. The control ligands were initially
docked with the binding sites of the five receptors, and the resulting
interactions were compared with standard reference ligands.

To prioritize new ligands for further studies, such as molecular dynamics and
others, we utilized various rules such as Lipinski, Pfizer, GSK, and golden
triangle rules. Additionally, important parameters for drug production,
including QED, SAscore, and MCE-18, were calculated using ADMETlab 2.024.

Table 3: The evaluation of GA-SVM-SVM:GA-SVM-SVM, GA-RF-RF:GA-QDA-SVM, GA-SVM-SVM:GA-QDA-SVM, and GA-RF-RF:GA-SVM-SVM models in the pipeline. | GA-SVM-SVM:
---
GA-SVM-SVM
precision | recall | f1-score | support
N/A | 0.75 | 0.74 | 0.74 | 461
EGFR+HER2 | 0.68 | 0.74 | 0.71 | 99
ER | 0.83 | 0.83 | 0.83 | 115
NF-$\kappa$B | 0.61 | 0.61 | 0.61 | 123
PR | 0.78 | 0.79 | 0.78 | 124
accuracy | 0.74 | 922
macro avg | 0.73 | 0.74 | 0.73 | 922
weighted avg | 0.74 | 0.74 | 0.74 | 922
| GA-RF-RF:
---
GA-QDA-SVM
precision | recall | f1-score | support
N/A | 0.73 | 0.76 | 0.74 | 461
EGFR+HER2 | 0.68 | 0.70 | 0.69 | 99
ER | 0.82 | 0.83 | 0.83 | 115
NF-$\kappa$B | 0.56 | 0.50 | 0.53 | 123
PR | 0.80 | 0.76 | 0.78 | 124
accuracy | 0.73 | 922
macro avg | 0.72 | 0.71 | 0.71 | 922
weighted avg | 0.72 | 0.73 | 0.73 | 922
| GA-SVM-SVM:
---
GA-QDA-SVM
precision | recall | f1-score | support
N/A | 0.75 | 0.74 | 0.74 | 461
EGFR+HER2 | 0.68 | 0.73 | 0.70 | 99
ER | 0.82 | 0.84 | 0.83 | 115
NF-$\kappa$B | 0.60 | 0.59 | 0.60 | 123
PR | 0.78 | 0.78 | 0.78 | 124
accuracy | 0.74 | 922
macro avg | 0.73 | 0.74 | 0.73 | 922
weighted avg | 0.74 | 0.74 | 0.74 | 922
| GA-RF-RF:
---
GA-SVM-SVM
precision | recall | f1-score | support
N/A | 0.73 | 0.76 | 0.74 | 461
EGFR+HER2 | 0.69 | 0.71 | 0.70 | 99
ER | 0.82 | 0.82 | 0.82 | 115
NF-$\kappa$B | 0.57 | 0.51 | 0.54 | 123
PR | 0.80 | 0.77 | 0.78 | 124
accuracy | 0.73 | 922
macro avg | 0.72 | 0.71 | 0.72 | 922
weighted avg | 0.73 | 0.73 | 0.73 | 922

##  3 Results and discussion

EGFR and HER2 receptors, shows 83.71% similarity in their residues using
sequence alignment with BLOSUM weight matrix and have a large similarity in
their 3D structure using Needleman-Wunsch alignment algorithm with BLOSUM-62
similarity matrix (Fig. S1 and S2). Moreover, near 70% of ligands in BindingDB
database with EGFR and HER2 targets, were identical. According to these
reasons, we merged two classes of EGFR and HER2 to just one EGFR/HER2 class.

The active/inactive and target classifiers utilized 128 and 64 selected
features, respectively (see tables S1 and S2). Various models were created
using the aforementioned methods. The top two models for the active/inactive
classifier were GA-SVM-SVM and GA-RF-RF, while for the target classifier, they
were GA-SVM-SVM and GA-QDA-SVM. The GA-SVM-SVM binary classifier was
constructed using the radial basis function (RBF) kernel with a gamma value of
0.1 and a regularization term of 1. Similarly, GA-RF-RF was built using the
Gini function with 400 trees and maximum depth until all leaves were pure. The
forest construction involved the use of bootstrap sampling. Furthermore, GA-
SVM-SVM and GA-QDA-SVM multi-classifiers were employed with the RBF kernel, a
gamma value of 0.01, and a regularization term of 10. The classification was
performed using the one-vs-one strategy, which has been shown to provide
higher prediction accuracy compared to the one-vs-rest approach25. As shown in
Table 2, the active/inactive classifiers GA-SVM-SVM and GA-RF-RF achieved
precision, recall, and F1-scores all above 0.7. Likewise, the target
classifiers GA-SVM-SVM and GA-RF-RF achieved precision, recall, and F1-scores
all above 0.9.

Table 4: Number of selected BindingDB molecules for each targets according to the threshold of 0, 80, 85, and 90% decision certainity for active/inactive prediction and the threshold of 0 and 90% decision certainity for target prediction. | Threshold
---|---
classes | 0:0 | 0:90 | 80:90 | 85:90 | 90:90
EGFR+HER2 | 172498 | 95123 | 19796 | 11068 | 4454
ER | 54101 | 22876 | 3613 | 2029 | 803
NF-$\kappa$B | 45452 | 16400 | 2499 | 1257 | 438
PR | 67323 | 14109 | 2300 | 1116 | 378
Figure 2: The ROC plots for different classes with one-vs-rest strategy.
Figure 3: The left plot shows the importance of features using permutation
importance method using GA-SVM-SVM model for target prediction and the right
one demonstrates hierarchical clustering dendrogram using pearson method to
find the correlation distance of each features. Figure 4: A simple
questionnaire dendrogram to separate ligands with a number of features and
determine the targets of them.

Subsequently, we generated a pipeline by combining different permutations of
the selected model. The precision, recall, f1-score, and support of these
models are presented in Table 3. Among the options, the GA-SVM-SVM:GA-SVM-SVM
model emerged as the most suitable pipeline, displaying superior performance
compared to others. This approach achieved an accuracy of 0.74 and an AUC of
0.94. Figure 2 showcases the ROC plots for each class using the one-vs-rest
strategy, further validating the effectiveness of the GA-SVM-SVM:GA-SVM-SVM
model for virtual screening. Table 4 provides insights into the number of
selected molecules from the BindingDB database for each target, based on
different predetermined thresholds. Notably, this table reveals the presence
of 4454, 803, 438, and 378 new inhibitor molecules for EGFR+HER2, ER,
NF-$\kappa$B, and PR, respectively. These novel inhibitors were selected with
90% precision in both the active/inactive and therapeutic classification
decision-making processes.

Figure 5: Docking results of new ligands obtained from virtual screening. The
pale dots in the following plots represent the active molecules in the
BindingDB database for each class, the filled dots represent the molecules
that participated in the construction of the model, and the red dots are the
new molecules proposed by the model obtained from the screening. Figure 6:
List of the molecules which are met all criteria.

In order to easily identify the target for each molecule, our objective was to
extract a straightforward rule from the model. To achieve this, we employed
permutation importance to determine the significance of each feature in the
model. Additionally, we utilized the Pearson method to create a hierarchical
clustering dendrogram, which helped us identify the correlation distance
(Euclidean distance) between features. Figure 3 displays the feature
importance and hierarchical clustering dendrogram. Based on these findings, we
constructed a simple questionnaire dendrogram for determining the target of
each molecule which are selected with 90% precision in both the
active/inactive and target classification decision-making, as illustrated in
Figure 4. The data presented in Figure 4 provide concise and effective
structure-activity relationship (SAR) information regarding the inhibitors.
For instance, NF-$\kappa$B inhibitors exhibit significantly lower values for
molecular walk count of order 10 (MWC10) and signal 10/weighted by
polarizability (Mor10p) compared to other inhibitors. Additionally, EGFR and
HER2 inhibitors demonstrate higher values of frequency of C–N at topological
distance 3 (F03[C-N]) and SHED Donor-Donor (SHED_DD) in comparison to ER and
PR inhibitors. The ratio of eigenvalue n.1 from augmented edge adjacency mat.
weighted by dipole moment (Eig01_AEA(dm)), radial distribution
function–105/weighted by polarizability (RDF105p), and radial distribution
function–105/weighted by I-state (RDF105s) differentiates PR and ER
inhibitors. These SAR information types effectively filter a significant
portion of large databases, thus accelerating early-stage drug discovery
projects that begin with extensive databases like GDB-1326. The classification
of molecules based on their therapeutic targets has garnered considerable
attention in the field of chemoinformatics27. These types of classifiers
expand on the concept of “Chemography”28, 29, which refers to the art of
navigating through chemical space. As evident from these figures, the
inhibitors cluster within specific regions of the selected chemical space,
aligning with the core objective of chemography.

In order to assess the binding energy of the molecules selected using the GA-
SVM-SVM:GA-SVM-SVM model, we employed molecular docking for both the chosen
molecules (with a multiplication of precision product exceeding 0.9 for both
active/inactive and target classification) and the molecule sets within each
class. The distribution of binding energy for these molecules, based on their
molecular weights, is depicted in Figure 5. In these plots, the pale dots
represent the active inhibitors labeled by the bindingDB database, while the
filled dots represent the active molecules utilized in constructing the GA-
SVM-SVM:GA-SVM-SVM model. The red dots correspond to new inhibitors, which
exhibit binding energy within the range of -15 to -5
$\mathrm{kcal}\text{/}\mathrm{mol}$. This range of binding energy proves to be
sufficiently suitable for forming protein-ligand complexes.

In order to prioritize further study on the new molecules, we applied several
medicinal criteria. The Lipinski rule suggests that ligands with a molecular
weight of less than or equal to 500 ($M_{W}\leq 500$), a logarithm of the
n-octanol/water distribution coefficient of less than or equal to 5 ($logP\leq
5$), a number of hydrogen bond acceptors of less than or equal to 10
($H_{acc}\leq 10$), and a number of hydrogen bond donors of less than or equal
to 5 ($H_{don}\leq 5$) exhibit good absorption or permeability. Accordingly,
376, 59, 91, and 35 ligands were accepted based on the Lipinski rule for
EGFR+HER2, ER, NF-$\kappa$B, and PR targets, respectively. The Pfizer rule
states that ligands with a $logP\geq 3$ and a topological polar surface area
of less than or equal to 75 ($TPSA\leq 75$) are likely to be toxic.
Consequently, 577, 30, 38, and 15 ligands passed the Pfizer rule for
EGFR+HER2, ER, NF-$\kappa$B, and PR targets, respectively. The GSK rule
suggests that ligands with a $MW\leq 400$ and $logP\leq 4$ may have a more
favorable ADMET (absorption, distribution, metabolism, excretion, and
toxicity) profile. Thus, 6, 11, 65, and 12 ligands were accepted based on the
GSK rule for EGFR+HER2, ER, NF-$\kappa$B, and PR targets, respectively.
Additionally, the golden triangle hypothesis proposes that ligands with a
$200\leq MW\leq 500$ and a logD (logarithm of the n-octanol/water distribution
coefficient at pH = 7.4) ranging from -2 to 5 ($-2\leq logD\leq 5$) may have a
more favorable ADMET profile. Consequently, 166, 53, 91, and 34 ligands
fulfilled the golden triangle criteria for EGFR+HER2, ER, NF-$\kappa$B, and PR
targets, respectively. Moreover, several parameters such as QED (desirability
functions based on eight drug-likeness related properties including $MW$,
$logP$, $N_{HBA}$, $N_{HBD}$, $PSA$, $N_{rotb}$, $N_{Ar}$), SAscore (synthetic
accessibility score based on a combination of fragment contributions and a
complexity penalty), and MCE-18 (medicinal chemistry evolution in 2018 score
molecules by novelty in terms of their cumulative sp3 complexity) were
considered favorable in the medical industry. Ligands with QED greater than
0.67, SAscore less than 6, and MCE-18 larger than 45 were deemed desirable.
Accordingly, 6, 14, 4, and 14 ligands met these criteria for EGFR+HER2, ER,
NF-$\kappa$B, and PR targets, respectively. The distributions of selected
molecules according to these rules are illustrated in Figures S3-S7. The
molecules depicted in Figure 6 satisfy all these criteria. Specifically, 2, 3,
3, and 8 ligands met all the parameters for EGFR+HER2, ER, NF-$\kappa$B, and
PR targets, respectively. As observed in this figure, each ligand within each
class exhibits unique structural properties.

##  4 Conclusion

In this research, we utilized various machine learning methods and employed
the GA-SVM-SVM and GA-RF-RF models for the active/inactive classification, as
well as GA-SVM-SVM and GA-QDA-SVM for therapeutic classification. Based on the
revenue generated by these methods in the pipeline, we selected GA-SVM-SVM:GA-
SVM-SVM as the most effective pipeline for virtual screening. This model
screened the BindingDB database inhibitors with varying precision ratios. As a
result, we identified 4454, 803, 438, and 378 new inhibitor molecules for
EGFR+HER2, ER, NF-$\kappa$B, and PR, respectively. Moreover, we introduced a
simple dendrogram to determine the target of each inhibitor with new ligands
which are predicted with over 90% precision in both active/inactive and target
prediction.

The molecules that exhibited a precision product exceeding 0.9 for both
active/inactive and target classification were selected for further evaluation
of their binding energy. Molecular docking analysis revealed that the binding
energies of these inhibitors to their therapeutic targets ranged from -15 to
-5 $\mathrm{kcal}\text{/}\mathrm{mol}$, which is considered suitable for
inhibiting the targets.

To prioritize further investigation of these new molecules, we applied several
filters, including the Lipinski, Pfizer, GSK, golden triangle rules, as well
as QED, SAscore, and MCE-18 parameters. Among these filters, 2, 3, 3, and 8
ligands met all the specified criteria for EGFR+HER2, ER, NF-$\kappa$B, and PR
targets, respectively. We believe that this study can provide valuable
insights for researchers working on the discovery of new inhibitors for breast
cancer.

The authors acknowledge the UNESCO UNISA iThemba-LABS/NRF Africa Chair in
Nanoscience & Nanotechnology (U2ACN2) and the Centre for High-Performance
Computing (CHPC), South Africa for providing computational resources and
facilities for this research project. The authors would like to express their
gratitude to all the members of the bioinformatics lab at Tarbiat Modares
University (TMU) for their valuable contributions in the form of discussions
and critical feedback on the manuscript.

Figures of sequence and structural alignment of EGFR and HER2 receptors,
distribution of ligands according to the Lipinski, Pfizer, GSK, golden
triangle rules and medical synthesis parameters and table of selected features
for active/inactive and target are presented in supporting information.

##  5 Data Availability

Scripts of this study is available on Github, allowing users to access the
source code. Visit our GitHub repository at https://github.com/parham-
rezaee/bindingdb-bc-virtual-screening.

## References

  * IsmaiL et al. 2018 IsmaiL, S.; UzaiRu, A.; Sagagi, B.; SuleiMan, M. S. Insilico Molecular Docking and Pharmacokinetic Studies of Selected Phytochemicals With Estrogen and Progesterone Receptors as Anticancer Agent for Breast Cancer. _Journal of the Turkish Chemical Society Section A: Chemistry_ 2018, _5_ , 1337–1350
  * Zhou et al. 2014 Zhou, Z.; Qiao, J. X.; Shetty, A.; Wu, G.; Huang, Y.; Davidson, N. E.; Wan, Y. RETRACTED ARTICLE: Regulation of estrogen receptor signaling in breast carcinogenesis and breast cancer therapy. _Cellular and Molecular Life Sciences_ 2014, _71_ , 1549–1549
  * Gnanaselvan et al. 2024 Gnanaselvan, S.; Yadav, S. A.; Manoharan, S. P. Structure-based virtual screening of anti-breast cancer compounds from Artemisia absinthium —insights through molecular docking, pharmacokinetics, and molecular dynamic simulations. _Journal of Biomolecular Structure and Dynamics_ 2024, _42_ , 3267–3285
  * Schwartz and Murray 2011 Schwartz, P. A.; Murray, B. W. Protein kinase biochemistry and drug discovery. _Bioorganic Chemistry_ 2011, _39_ , 192–210
  * Yang et al. 2015 Yang, J. C.-H. et al. Afatinib versus cisplatin-based chemotherapy for EGFR mutation-positive lung adenocarcinoma (LUX-Lung 3 and LUX-Lung 6): analysis of overall survival data from two randomised, phase 3 trials. _The Lancet Oncology_ 2015, _16_ , 141–151
  * Johnston, S.R.D. and Leary, A. 2006 Johnston, S.R.D.; Leary, A. Lapatinib: A novel EGFR/HER2 tyrosine kinase inhibitor for cancer. _Drugs of Today_ 2006, _42_ , 441
  * Brennan et al. 2000 Brennan, P. J.; Kumogai, T.; Berezov, A.; Murali, R.; Greene, M. I. HER2/Neu: mechanisms of dimerization/oligomerization. _Oncogene_ 2000, _19_ , 6093–6101
  * Li et al. 2014 Li, L.; Wang, Y.; Qi, B.; Yuan, D.; Dong, S.; Guo, D.; Zhang, C.; Yu, M. Suppression of PMA-induced tumor cell invasion and migration by ginsenoside Rg1 via the inhibition of NF-κB-dependent MMP-9 expression. _Oncology Reports_ 2014, _32_ , 1779–1786
  * Purawarga Matada et al. 2022 Purawarga Matada, G. S.; Dhiwar, P. S.; Abbas, N.; Singh, E.; Ghara, A.; Das, A.; Bhargava, S. V. Molecular docking and molecular dynamic studies: screening of phytochemicals against EGFR, HER2, estrogen and NF-KB receptors for their potential use in breast cancer. _Journal of Biomolecular Structure and Dynamics_ 2022, _40_ , 6183–6192
  * Asadollahi-Baboli and Mani-Varnosfaderani 2013 Asadollahi-Baboli, M.; Mani-Varnosfaderani, A. Shuffling multivariate adaptive regression splines as a predictive method for modeling of novel pyridylmethylthio derivatives as VEGFR2 inhibitors. _Medicinal Chemistry Research_ 2013, _22_ , 2645–2653
  * Zubair and Frieri 2013 Zubair, A.; Frieri, M. Role of Nuclear Factor-$\kappa$B in Breast and Colorectal Cancer. _Current Allergy and Asthma Reports_ 2013, _13_ , 44–49
  * Alberti et al. 2012 Alberti, C.; Pinciroli, P.; Valeri, B.; Ferri, R.; Ditto, A.; Umezawa, K.; Sensi, M.; Canevari, S.; Tomassetti, A. Ligand-dependent EGFR activation induces the co-expression of IL-6 and PAI-1 via the NFkB pathway in advanced-stage epithelial ovarian cancer. _Oncogene_ 2012, _31_ , 4139–4149
  * Peralta et al. 2009 Peralta, E. A.; Murphy, L. L.; Minnis, J.; Louis, S.; Dunnington, G. L. American Ginseng Inhibits Induced COX-2 and NFKB Activation in Breast Cancer Cells. _Journal of Surgical Research_ 2009, _157_ , 261–267
  * Sun et al. 2022 Sun, D.; Gao, W.; Hu, H.; Zhou, S. Why 90% of clinical drug development fails and how to improve it? _Acta Pharmaceutica Sinica B_ 2022, _12_ , 3049–3062
  * Alamdari et al. 2012 Alamdari, R.; Mani-Varnosfaderani, A.; Asadollahi-Baboli, M.; Khalafi-Nezhad, A. Monte Carlo sampling and multivariate adaptive regression splines as tools for QSAR modelling of HIV-1 reverse transcriptase inhibitors. _SAR and QSAR in Environmental Research_ 2012, _23_ , 665–682
  * Bajorath 2015 Bajorath, J. Computer-aided drug discovery. _F1000Research_ 2015, _4_ , 630
  * Sadybekov and Katritch 2023 Sadybekov, A. V.; Katritch, V. Computational approaches streamlining drug discovery. _Nature_ 2023, _616_ , 673–685
  * Jalali-Heravi and Mani-Varnosfaderani 2011 Jalali-Heravi, M.; Mani-Varnosfaderani, A. QSAR modelling of integrin antagonists using enhanced Bayesian regularised genetic neural networks. _SAR and QSAR in Environmental Research_ 2011, _22_ , 293–314
  * Gilson et al. 2016 Gilson, M. K.; Liu, T.; Baitaluk, M.; Nicola, G.; Hwang, L.; Chong, J. BindingDB in 2015: A public database for medicinal chemistry, computational chemistry and systems pharmacology. _Nucleic Acids Research_ 2016, _44_ , D1045–D1053
  * O’Boyle et al. 2011 O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. _Journal of Cheminformatics_ 2011, _3_ , 33
  * Frisch 2009 Frisch, A. _Gaussian 09: user’s reference_ ; Gaussian: Wallingford, Conn., 2009; OCLC: 711965588
  * Mauri 2020 Mauri, A. In _Ecotoxicological QSARs_ ; Roy, K., Ed.; Springer US: New York, NY, 2020; pp 801–820
  * Trott and Olson 2010 Trott, O.; Olson, A. J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. _Journal of Computational Chemistry_ 2010, _31_ , 455–461
  * Xiong et al. 2021 Xiong, G.; Wu, Z.; Yi, J.; Fu, L.; Yang, Z.; Hsieh, C.; Yin, M.; Zeng, X.; Wu, C.; Lu, A.; Chen, X.; Hou, T.; Cao, D. ADMETlab 2.0: an integrated online platform for accurate and comprehensive predictions of ADMET properties. _Nucleic Acids Research_ 2021, _49_ , W5–W14
  * Jalali‐Heravi et al. 2013 Jalali‐Heravi, M.; Mani‐Varnosfaderani, A.; Valadkhani, A. Integrated One‐Against‐One Classifiers as Tools for Virtual Screening of Compound Databases: A Case Study with CNS Inhibitors. _Molecular Informatics_ 2013, _32_ , 742–753
  * Blum and Reymond 2009 Blum, L. C.; Reymond, J.-L. 970 Million Druglike Small Molecules for Virtual Screening in the Chemical Universe Database GDB-13. _Journal of the American Chemical Society_ 2009, _131_ , 8732–8733
  * Jalali-Heravi et al. 2011 Jalali-Heravi, M.; Mani-Varnosfaderani, A.; Jahromi, P. E.; Mahmoodi, M. M.; Taherinia, D. Classification of anti-HIV compounds using counterpropagation artificial neural networks and decision trees. _SAR and QSAR in Environmental Research_ 2011, _22_ , 639–660
  * Von Korff and Hilpert 2006 Von Korff, M.; Hilpert, K. Assessing the Predictive Power of Unsupervised Visualization Techniques to Improve the Identification of GPCR-Focused Compound Libraries. _Journal of Chemical Information and Modeling_ 2006, _46_ , 1580–1587
  * Jalali‐Heravi and Mani‐Varnosfaderani 2012 Jalali‐Heravi, M.; Mani‐Varnosfaderani, A. Navigating Drug‐Like Chemical Space of Anticancer Molecules Using Genetic Algorithms and Counterpropagation Artificial Neural Networks. _Molecular Informatics_ 2012, _31_ , 63–74


# Differential equation for the flow rate of discharging silos based on energy
balance

J. R. Darias1, Marcos A. Madrid2,3, Luis A. Pugnaloni4 1Laboratorio de Óptica
y Fluidos, Universidad Simón Bolívar, Apartado Postal 89000, Caracas 1080-A,
Venezuela.
2Departamento de Ingeniería Mecánica, Facultad Regional La Plata, Universidad
Tecnológica Nacional, CONICET, Av. 60 Esq. 124, 1900 La Plata, Argentina.
3Instituto de Física de Líquidos y Sistemas Biológicos (CONICET La Plata,
UNLP), Calle 59 Nro 789, 1900 La Plata, Argentina.
4Departamento de Física, Facultad de Ciencias Exactas y Naturales, Universidad
Nacional de La Pampa, CONICET, Uruguay 151, 6300 Santa Rosa (La Pampa),
Argentina.

###### Abstract

Since the early work of Hagen in 1852 and Beverloo et al. in 1961, the flow
rate of granular material discharging through a circular orifice from a silo
has been described by means of dimensional analysis and experimental fits, and
explained through the “free fall arch” model. Here, in contrast with the
traditional approach, we derive a differential equation based on the energy
balance of the system. This equation is consistent with the well known
Beverloo rule thanks to a compensation of energy terms. Moreover, this new
equation can be used to explore new conditions for silo discharges. In
particular, we show how the effect of friction on the flow rate can be
predicted. The theory is validated using discrete element method simulations.

###### pacs:

45.70.-n, 45.70.Mg

##  I Introduction

When a granular material (such as seeds) flow through an orifice at the base
of a silo, the resulting flow rate presents peculiar features when compared
with the better known phenomenology of inviscid fluids. In particular, the
flow rate of grains does not depend on the height $h$ of the material in the
silo. This was early noted by Hagen in 1852 hagen1852  (for a translation of
the German publication see Ref. tighe2007 ). Moreover, Hagen showed that the
mass flow rate $Q$ scales as $D_{\textrm{o}}^{5/2}$, with $D_{\textrm{o}}$ the
diameter of the circular orifice. Compare this with $D_{\textrm{o}}^{2}$,
which is valid for inviscid fluids. Hagen also provided the first heuristic
postulate to explain this 5/2 power which was later termed the “free fall
arch” model BrownBook . In brief, Hagen postulated that grains move downwards
in the silo at a very low velocity until they arrive at a region (about one
orifice radius tall) close to the opening at which the local density is low
and grains perform a simple free fall. From here, the typical vertical
velocity $v$ at which grains cross the plane of the orifice can be estimated
($v_{\textrm{out}}=\sqrt{gD_{\textrm{o}}}$, with $g$ the acceleration of
gravity) and the mas flow rate calculated as
$Q=\pi(D_{\textrm{o}}/2)^{2}\rho_{\textrm{o}}v_{\textrm{out}}=(\pi/4)\rho_{\textrm{o}}\sqrt{g}D_{\textrm{o}}^{5/2}$,
being $\rho_{\textrm{o}}$ the local apparent density at the orifice. It is
customary to replace $\rho_{\textrm{o}}$ by the density in the bulk of the
silo $\rho_{\textrm{b}}$. Therefore, the $\pi/4$ coefficient is replaced by a
constant $C$ that is later fitted to the experimental data. Hagen also noted
that this expression does not agree entirely well with the data. The solution
to this was to consider that the effective orifice is about one grain
diameter, $d$, smaller due to boundary effects which leads to
$Q=C\rho_{b}\sqrt{g}(D_{\textrm{o}}-d)^{5/2}$.

In 1961 Beverloo et al. carried out a series of experiments and found a more
suitable expression to fit the mass flow rate beverloo . This expression (see
below) became widely accepted, and is usually named as the Beverloo rule.
Beverloo et al. did not refer to Hagen’s work, but found the 5/2 power based
on dimensional analysis. Others in the decade preceding Beverloo proposed
expression that were less successful. The improvement with respect to Hagen
was the introduction of an additional constant $k$ to the boundary effect
correction. The Beverloo rule states that the mass flow rate is

$Q=C\rho_{\textrm{b}}\sqrt{g}(D_{\textrm{o}}-k\,d)^{5/2},$ (1)

where $k$ and $C$ are two fitting dimensionless constants. The fitted value
for $k_{\textrm{experim}}$ may vary between $1.4$ and $3$ depending on the
shape and size dispersion of the grains. However, $C_{\textrm{experim}}\approx
0.58$ with almost no influence of the type of material the grains are made of
NeddermanBook . The origin for this “universal” value of $C$ did not receive
much attention until recently madrid2018 .

The problem of discharge of grains through an orifice has been revisited by a
number of authors (see for example nedderman1982 ; tuzun ; savage ; ristow ;
mankoc2007 ; kondic2014 ; wilson2014 ; goldberg2016 ; villagran  and
references therein). The basic phenomenology has been confirmed in all
studies, i.e.: (i) the flow rate does not depend on the column height, (ii)
the flow rate scales with $D_{\textrm{o}}^{5/2}$, (iii) the prefactor is
$C\approx 0.58$ for almost all materials. However, Rubio-Largo et al. provided
evidence that the concept of “free fall arch” may not be a realistic picture
of the internal dynamics in the silo around the orifice rubio2015 .

This mechanical problem of discharge should be described from first principles
via energy balance as it is done for a fluid (e.g., Bernoulli’s law). To our
knowledge, this has not been done successfully so far, despite some attempts
(see for example mcdougall1965 ). One possible reason for this is that the
rheological response of the grains while flowing in a silo has received less
attention, leaving a gap in a basic component of the energy balance analysis:
_the dissipated power_. In this direction, Staron et al. staron2012  have
proved that the flow in a discharging silo can be described by introducing the
$\mu(I)$-rheology model dacruz2005  into the Navier–Stokes equations. This is,
in practice, the introduction of an effective dissipation term into the
equations. Recently, we have shown that, at a global scale in the silo, the
flow is consistent with the quasistatic limit in the $\mu(I)$-rheology
madrid2018 . In that work, we have provided an expression for the energy
dissipated during the discharge of a silo that will be revisited in what
follows.

In this work, we derive a differential equation for $Q$ as a function of time
during the discharge. This is achieved by using the work–energy theorem for
the system of interest, which is defined as the set of grains that remain
inside the silo at any given time. To calculate the energy dissipation we use
a revised expression to the one presented in Ref. madrid2018 . The result is
consistent with the Beverloo rule. We provide a theoretical estimate for the
value of $C$ that is remarkably close to the value obtained by experimental
fits. Moreover, the new equation predicts an increase in the flow rate when
the friction coefficient is reduced, which is consistent with simulation
results.

##  II System definitions

Figure 1: Sketch of the axial cross section of a cylindrical silo.

We consider a cylindrical silo (see Fig. 1) of diameter $D_{\rm s}$ (radius
$R_{\textrm{s}}=D_{\rm s}/2$ and cross section $A_{\textrm{s}}=\pi
R_{\textrm{s}}^{2}$) with a flat base. This base has an orifice of diameter
$D_{\textrm{o}}$ in its center (radius $R_{\textrm{o}}=D_{\textrm{o}}/2$ and
cross section $A_{\textrm{o}}=\pi R_{\textrm{o}}^{2}$). With the orifice
covered, the silo is filled with an initial mass $M_{\textrm{ini}}$ of a
granular material that fills the silo up to a certain height. We will consider
that the bulk apparent density of the material $\rho_{\textrm{b}}$ is
homogeneous throughout the column as a first order approximation. The mass of
each grain is denoted by $m$.

When the orifice is opened, the discharge of grains starts. During the
discharge, the mass $M(t)$ inside the silo at time $t$ can be written as

${}M(t)=\rho_{\textrm{b}}A_{\textrm{s}}z(t)=2\rho_{\textrm{b}}A_{\textrm{s}}z_{\textrm{cm}}(t),$
(2)

where $z(t)$ is the height of the column of grains, and
$z_{\textrm{cm}}(t)=z(t)/2$ is the center of mass of the granular column at
time $t$. Therefore, the mass flow rate $Q(t)$ is

$Q(t)=-\dot{M}(t)=2\rho_{\textrm{b}}A_{\textrm{s}}v_{\textrm{cm}}(t),$ (3)

where $v_{\textrm{cm}}(t)=|\dot{z}_{\textrm{cm}}(t)|$ is the speed of the
center of mass of the granular column. Note that since $M(t)$ decreases with
time during the discharge, $\dot{M}$ and $\dot{z}_{\textrm{cm}}$ are negative.

##  III Simulations

To validate our theoretical model we carried out a series of DEM simulations
of spherical particles in a cylindrical silo as described in the previous
section. We used the LIGGGHTS liggghts  implementation with a
particle–particle Hertz interaction and Coulomb criterion using a Young
modulus $Y=70$ MPa, Poison ratio $\nu=0.25$, restitution coefficient $e=0.95$
and friction coefficient $0.1<\mu<1.0$. The same interaction applies for the
particle–walls contacts. Details on the particle–particle interactions are
given in Appendix A. The particle diameters are $d=1$ mm and their material
density is $\rho=2500\ \rm{kg/m}^{3}$. The silo diameter is
$D_{\mathrm{s}}=30\ \rm{mm}$ (some test have been run also with
$D_{\mathrm{s}}=24\ \rm{mm}$). The orifice diameter is varied in the range
$6.0\ {\rm mm}\leq D_{\rm o}\leq 10.0\ {\rm mm}$. Particles are poured into
the silo to fill a height $z(t=0)\approx 10D_{\mathrm{s}}$ (which implies up
to $3\times 10^{5}$ grains, depending on the silo diameter). The orifice is
initially blocked by a plug. We let the grains come to rest in the silo by
waiting until the kinetic energy per particle falls below $10^{-10}$ J. Then,
we remove the plug and allow the material to discharge through the orifice.
Although particles are monosized, we do not observe crystalline structures in
the simulations. The magnitude of the acceleration of gravity is set to
$g=9.81$ m/s2 and the integration time step is $\Delta t=5\times 10^{-6}$ s.

##  IV Work–energy theorem and dissipated power

We focus on the system composed of the grains inside the silo at any given
time. According to the work–energy theorem, the change in kinetic energy
$\dot{K}_{\textrm{in}}$ of this system of grains is

$\dot{K}_{\textrm{in}}=\dot{W}_{\textrm{g}}-\dot{W}_{\textrm{out}}+\dot{W}_{\textrm{el}}-\dot{W}_{\textrm{d}},$
(4)

where $\dot{W}_{\textrm{g}}$ is the power injected by the force of gravity
acting on the grains, $\dot{W}_{out}$ is the power loss due to the grains that
leave the silo through the orifice at a velocity $v_{\textrm{out}}$,
$\dot{W}_{\textrm{el}}$ is the “elastic power”, i.e., the rate of change of
the elastic energy of the grain–grain contacts, and $\dot{W}_{\textrm{d}}$ is
the dissipated power due to the non-conservative interactions (friction and
inelastic collisions between the grains and between the grains and the walls).
Some of the terms in Eq. (4) are in fact negligible ($\dot{K}_{\textrm{in}}$
and $\dot{W}_{\textrm{el}}$) and some are easy to calculate from basic
mechanics ($\dot{W}_{\textrm{g}}$ and $\dot{W}_{\textrm{out}}$). We discuss
those contributions in Appendix B. In the remaining of this subsection we
focus on the dissipated power $\dot{W}_{\textrm{d}}$.

In a recent work madrid2018 , we have shown that the power dissipated during a
silo discharge can be calculated by assuming that the flow is consistent with
a quasistatic shear flow in the framework of the $\mu(I)$-rheology dacruz2005
. Here, we present an improved expression for $\dot{W}_{\textrm{d}}$ in which
we consider in more detail the local dissipation due to the convergent flow
region in the bottom part of the silo.

Figure 2: Sketch of a simple plane shear cell. The top plate is driven at
constant velocity $v$.

Let us first consider a simple plane shear configuration (see Fig. 2).
According to the $\mu(I)$-rheology model dacruz2005 , the average tangential
stress $\sigma_{xz}$ required to keep the top plate moving at constant
velocity $v$ can be written as $\sigma_{xz}=\mu(I)\sigma_{zz}$. Here,
$\sigma_{zz}$ is the confining stress and $\mu(I)$ is an “effective friction
coefficient” that accounts for all the complex interactions in the granular
sample between the plates. The effective friction depends on the inertial
number $I=vd\sqrt{\rho}/(L\sqrt{\sigma_{zz}})$, with $\rho$ the density of the
material the grains are made of dacruz2005 . The inertial number suffices to
characterize the flow as long as the grains are stiff and $L\gg d$.

The power $\dot{W}_{\rm d}$ dissipated during the motion of the top plate in
Fig. 2 is simply given by $\dot{W}_{\rm d}=\sigma_{xz}Av$, where $A$ is the
total area of the top plate. Hence, the dissipated power can be written as
$\dot{W}_{\rm d}=\mu(I)\sigma_{zz}Av$. In this expression, the effect of the
properties of the granular material on the dissipation comes only through
$\mu(I)$ since $\sigma_{zz}$, $A$ and $v$ are control variables of the
experiment. Interestingly, it has been shown that $\mu(I)$ follows a universal
curve for all granular materials as long as the particle–particle friction is
high enough (roughly above $0.4$) dacruz2005 . However, this curve depends on
geometrical factors like dimensionality azema2014 . In the quasistatic limit
(i.e., $I\ll 1$, with $I<0.01$ being a typical criterion used for practical
proposes), the value of $\mu(I=0)$ becomes unique to all granular materials
for a given geometry if the particle–particle friction is above $0.4$.
Therefore, the dissipated power becomes independent from the details of the
particle–particle interactions in this limit (i.e., when $I\rightarrow 0$ and
$\mu>0.4$).

We postulate that the flow inside a discharging silo can be modeled in a
similar way as in the simple plane shear geometry. In the cylindrical silo
geometry, up and away from the converging flow observed at the bottom of the
silo, the confining pressure can be taken as the mean radial pressure
$\langle\sigma_{rr}\rangle(t)$ (where the angular brackets indicate an average
over the height of the granular column), the driving velocity as the velocity
of the free surface of the granular column $v=2v_{\textrm{cm}}$ and the “plate
area” as the area of contact between the grains and the silo lateral walls
[$A_{1}(t)=\pi D_{\rm s}2z_{\textrm{cm}}(t)$]. Note that $v$ is only a
velocity that is characteristic for the motion of the system at the
macroscopic scale (as the plate velocity in the plane shear experiment) and
does not need to be compared with the actual velocities of the grains nor with
the velocity gradients in the system.

Close to the bottom of the silo, over a height of about $R_{\rm s}$, the flow
converges to the orifice and the relevant confining pressure is no longer
radial in this section of the silo. We therefore use, as a proxy for the
confining pressure in this region, the vertical component $\sigma^{\rm
bott}_{zz}(t)$ averaged over the bottom part of the silo (from the base to a
height $R_{\rm s}$). The lateral area of this region is $A=\pi D_{\rm s}R_{\rm
s}$, and we use also the characteristic velocity $v=2v_{\textrm{cm}}$.
However, there exist a small region of height $\approx R_{\textrm{o}}$ right
above the base that does not contribute to the dissipation of energy. This is
due to the fact that grains close to the solid base do not move significantly
and then do not dissipate energy. Of course, grains in the neighborhood of the
orifice do move at high velocities, but here the packing fraction is so low
that there are very few particle–particle interactions to dissipate energy.
Therefore, for the converging flow zone we take the effective lateral area as
$A_{2}=2\pi R_{\rm s}(R_{\rm s}-D_{\rm o}/2)$. Using these approximations we
can write the dissipated power in the entire granular column as

$\dot{W}_{\textrm{d}}(t)=\mu(I)[\langle\sigma_{rr}\rangle(t)A_{1}(t)+\alpha\sigma^{\rm
bott}_{zz}(t)A_{2}]2v_{\textrm{cm}}(t).$ (5)

The first term in Eq. (5) accounts for the dissipation along the upper part of
the silo, whereas the second term accounts for the contribution of the
converging flow. The constant $\alpha$ in this second term is introduced as a
correction factor since the active area, characteristic velocity and confining
pressure in the converging zone are only estimates that should give the
correct scaling on the various variables but not necessarily the correct
numerical value for the dissipation.

Using Eqs. (2) and (3), we can write Eq. (5) as

$\displaystyle\dot{W}_{\textrm{d}}(t)=$ $\displaystyle-\mu(I)2\pi R_{\rm
s}\left[\langle\sigma_{rr}\rangle(t)\frac{M(t)}{\rho_{\rm b}A_{\rm s}}\right.$
$\displaystyle+\left.\alpha\sigma^{\rm bott}_{zz}(t)\left(R_{\rm
s}-\frac{D_{\rm o}}{2}\right)\right]\frac{\dot{M}(t)}{\rho_{\rm b}A_{\rm s}}.$
(6)

##  V A differential equation for the mass flow rate

Collecting all terms for the energy balance from Appendix B and the previous
section, i.e., plugging Eqs. (30), (32), (34), (35) and (6) into Eq. (4), one
can obtain a differential equation for the mass flow rate
$Q\equiv-\dot{M}(t)$. As we mentioned in Appendix B, some of these terms are
in fact negligible. In particular, we take $\dot{K}_{\textrm{in}}=0$ and
$\dot{E}_{\textrm{el}}=0$. Finally, Eq. (4) can be written as

$\displaystyle 0=$ $\displaystyle-\frac{g}{\rho_{\rm b}A_{\rm
s}}M(t)\dot{M}(t)+\frac{\dot{M}^{3}(t)}{2\rho_{\rm o}^{2}A_{\rm
o}^{2}}+\mu(I)2\pi R_{\rm s}$ (7)
$\displaystyle\times\left[\langle\sigma_{rr}\rangle(t)\frac{M(t)}{\rho_{\rm
b}A_{\rm s}}+\alpha\sigma^{\rm bott}_{zz}(t)\left(R_{\rm s}-\frac{D_{\rm
o}}{2}\right)\right]\frac{\dot{M}(t)}{\rho_{\rm b}A_{\rm s}}.$

Solving for $-\dot{M}$ foot

$\displaystyle-\dot{M}(t)=$ $\displaystyle\frac{\pi\sqrt{2}}{4}\rho_{\rm
o}\sqrt{g}D_{\rm o}^{2}\left\\{\frac{M(t)}{\rho_{\rm b}A_{\rm
s}}-\frac{\mu(I)2\pi R_{\rm s}}{g\rho_{\rm b}A_{\rm s}}\right.$
$\displaystyle\times$
$\displaystyle\left.\left[\langle\sigma_{rr}\rangle(t)\frac{M(t)}{\rho_{\rm
b}A_{\rm s}}+\alpha\sigma^{\rm bott}_{zz}(t)(R_{\rm s}-\frac{D_{\rm
o}}{2})\right]\right\\}^{1/2}.$

Equation (V) is a first order differential equation for $M(t)$ that can be
closed with an initial condition such as $M(t=0)=M_{0}$. To solve this
equation it is necessary to know $\langle\sigma_{rr}\rangle(t)$, $\sigma^{\rm
bott}_{zz}(t)$, $\mu(I)$ and $\alpha$. In the next section we revise the
Walters model for the pressure in discharging silos to obtain analytical
expressions for $\langle\sigma_{rr}\rangle(t)$ and $\sigma^{\rm bott}_{zz}(t)$
to close Eq. (V).

It is worth mentioning at this point that Eq. (7) reduces to the equation for
an inviscid fluid if the last term that accounts for the dissipated power is
neglected. In this case we obtain

$\displaystyle-\frac{\dot{M}^{2}(t)}{\rho^{2}_{\textrm{o}}A^{2}_{\textrm{o}}}$
$\displaystyle=$ $\displaystyle
2g\frac{M(t)}{\rho_{\textrm{b}}A_{\textrm{s}}},$ $\displaystyle v_{\rm out}$
$\displaystyle=$ $\displaystyle\sqrt{2gz(t)},$ (9)

where $z(t)$ is the column height and $v_{\rm
out}=\dot{M}(t)/(\rho^{2}_{\textrm{o}}A^{2}_{\textrm{o}})$ is the mean
velocity of the outflowing material.

##  VI Pressure in a discharging silo

###  VI.1 Walters dynamic stresses

Walters walters1973  and Walker walker1966  have developed models for the
pressure during silo discharge following an approach similar to Janssen’s. In
order to close Eq. (V) we will use these previous developments. The expression
for $\sigma_{zz}$ and $\sigma_{rr}$ at a given depth $z^{\prime}$ (measured
from the free surface of the granular column) into the moving column is

$\displaystyle\sigma_{zz}^{\rm Walters}(z^{\prime})$
$\displaystyle=\frac{g\rho_{\rm b}D_{\rm s}}{4B}[1-{\rm
e}^{-4Bz^{\prime}/D_{\rm s}}],$ (10) $\displaystyle\sigma_{rr}^{\rm
Walters}(z^{\prime})$ $\displaystyle=\frac{g\rho_{\rm b}D_{\rm
s}}{4\tan\phi}[1-{\rm e}^{-4Bz^{\prime}/D_{\rm s}}],$ (11)

where $\tan\phi$ is the effective friction corresponding to the wall yield
locus NeddermanBook , and $B$ plays the role of the well-known Janssen’s force
redirection factor. For discharging (not static) silos Walters obtains
walters1973

$\displaystyle B$
$\displaystyle=\frac{\tan\phi\cos^{2}\delta}{(1+\sin^{2}\delta)-2y\sin\delta},$
(12) $\displaystyle y$ $\displaystyle=\frac{2}{3c}[1-(1-c)^{3/2}],$
$\displaystyle c$ $\displaystyle=\frac{\tan^{2}\phi}{\tan^{2}\delta}.$

Here, $\phi$ is the effective friction angle for the wall yield locus and
$\delta$ is the effective friction angle for the internal yield locus. As
discussed by Nedderman (see section 3.7 in Ref. NeddermanBook ), the wall
yield locus $\phi$ should not be set simply as ${\rm atan}(\mu_{\rm wall})$.
Instead, the Jenike’s rule should be applied

$\tan\phi=\left\\{{\begin{array}[]{l}\mu_{\rm wall}\text{ if
}\sin\delta\geq\mu_{\rm wall}\\\ \sin\delta\text{ if }\sin\delta<\mu_{\rm
wall}\end{array}}\right..$ (13)

Figure 3: (a) $\langle\sigma_{rr}\rangle$ as a function of the total mass $M$
in the silo. (b) $\sigma^{\rm bott}_{zz}$ as a function of the total mass $M$
in the silo. Symbols correspond to DEM simulations for different orifice
diameters $D_{\rm o}$ and different silo diameters $D_{\rm s}$. Unless
otherwise stated $D_{\rm s}=30d$. The particle–particle and particle–wall
friction coefficients were set to $\mu=0.5$. The solid lines correspond to:
(a) $\langle\sigma_{rr}\rangle$ according to Eq. (VI.1) and (b) $\sigma^{\rm
bott}_{zz}(z)$ given by Eq. (14). We have set $\tan\delta=0.204$ to fit
simultaneously both pressures. For $\tan\phi$ we used $\tan\phi=\sin\delta$,
following Eq. (13).

For the vertical pressure on the converging flow zone we will use as an
estimate

$\displaystyle\sigma^{\rm bott}_{zz}(z)$ $\displaystyle=\sigma_{zz}^{\rm
Walters}(z^{\prime}=z)$ (14) $\displaystyle=\frac{g\rho_{\rm b}D_{\rm
s}}{4B}[1-{\rm e}^{-4Bz/D_{\rm s}}],$

being $z$ the granular column height. Note that Walters defines $z^{\prime}$
in the negative direction, therefore $\sigma_{zz}^{\rm Walters}(z^{\prime}=z)$
is simply the vertical pressure on the base of the silo.

Since in Eq. (V) we also require the average $\langle\sigma_{rr}\rangle$ over
the entire column of grains, we average Eq. (11) from $z^{\prime}=0$ to the
total depth $z^{\prime}=z$ of the column, which yields

$\displaystyle\langle\sigma_{rr}\rangle(z)$
$\displaystyle=\frac{1}{z}\int_{0}^{z}{\sigma_{rr}^{\rm
Walters}(z^{\prime})dz^{\prime}}$ $\displaystyle=\frac{g\rho_{\rm b}D_{\rm
s}}{4\tan\phi}\left[1+\frac{{\rm e}^{-4Bz/D_{\rm s}}-1}{4Bz/D_{\rm
s}}\right].$ (15)

In Fig. 3 we show $\sigma^{\rm bott}_{zz}$ and $\langle\sigma_{rr}\rangle$
obtained from DEM simulations along with the theoretical predictions form Eq.
(14) and (VI.1), respectively. The particle–wall friction was set $\mu_{\rm
wall}=0.5$, which is larger than the internal friction of the material. The
curves have been fitted setting $\tan\delta=0.204$ and $\tan\phi=\sin\delta$,
following Eq. (13).

There are two interesting features to emphasize in Fig. 3. Firstly,
$\langle\sigma_{rr}\rangle$ saturates much more slowly than the bottom
pressure $\sigma^{\rm bott}_{zz}$. Secondly, the prediction for
$\langle\sigma_{rr}\rangle$ fails to some extent for low column heights.
Although this could be improved, we will show that these expression for the
pressure are sufficient to obtain new valuable insights into the silo
discharge.

###  VI.2 Asymptotic pressure

A first order approximation to the solution for Eq. (V) can be obtained by
replacing the asymptotic limit of Eqs. (14) and (VI.1) for high columns ($z\gg
D_{\rm s}$). These asymptotic expressions are

$\displaystyle\sigma^{\rm bott}_{zz}(z)$ $\displaystyle=\frac{g\rho_{\rm
b}D_{\rm s}}{4B}[1+O(z^{2})]$ (16) $\displaystyle\langle\sigma_{rr}\rangle(z)$
$\displaystyle=\frac{g\rho_{\rm b}D_{\rm s}}{4\tan\phi}\left[1-\frac{D_{\rm
s}}{4Bz}+O(z^{2})\right]$

##  VII Comparison with the Beverloo rule

Figure 4: Dissipated power $\dot{W}_{\rm d}$ (scaled by the mass flow rate
$Q$) as a function of the total mass $M$ in the silo. Symbols correspond to
DEM simulations for different orifice diameters $D_{\rm o}$ and different silo
diameters $D_{\rm s}$. If not stated $D_{\rm s}=30d$. The black solid line
corresponds to the expression proposed in Eq. (6). The blue dashed and orange
dot-dashed lines correspond, respectively, to the first and second term in Eq.
(6) using $\mu(I)=\tan\phi=\sin\delta=0.2$ and $\alpha=2.5$ as discussed in
the text.

###  VII.1 The asymptotic equation and the $5/2$ power law

Replacing Eq. (16) into Eq. (V), and using $z=M(t)/(\rho_{\rm b}A_{\rm s})$
(see Eq. (2)) we obtain

$\displaystyle-\dot{M}(t)=$ $\displaystyle\frac{\pi\sqrt{2}}{4}\rho_{\rm
o}\sqrt{g}D_{\rm o}^{2}\left[\frac{M(t)}{\rho_{\rm b}A_{\rm
s}}\left(1-\frac{\mu(I)}{\tan\phi}\right)\right.$ (17)
$\displaystyle+\left.\frac{\mu(I)R_{\rm
s}}{2B}\left(\frac{1}{\tan\phi}-2\alpha\right)+\frac{\alpha\mu(I)}{2B}D_{\rm
o}\right]^{1/2}$

Equation (17) becomes similar to the Hagen’s expression for the flow rate if
we select $\mu(I)=\tan\phi=(2\alpha)^{-1}$ and $B=1/4$. Under these conditions

$-\dot{M}(t)=\frac{\pi\sqrt{2}}{4}\rho_{\rm o}\sqrt{g}D_{\rm o}^{2.5}.$ (18)

As we showed in the previous section, the Walters expressions for the pressure
fit the DEM data if $\tan\delta=0.204$ and $\tan\phi=\sin\delta=0.20$. For
these values, Eq. (12) yields $B=0.26\approx 1/4$. Therefore, the well known
Hagen’s relation will hold if we simply set $\mu(I)=0.20$ and $\alpha=2.5$.

To validate the values of $\mu(I)$ and $\alpha$ required for Eq. (17) to
reduce to Eq. (18), we have calculated from our DEM simulations the energy
dissipated during discharge (see Appendix A for details). In Fig. 4 we plot
the dissipated power along with the expression provided by Eq. (6) using
$\mu(I)=0.20$ and $\alpha=2.5$. As we can see, the agreement is remarkable not
only for high columns (where the asymptotic limit was used to predict $\mu(I)$
and $\alpha$) but also for low column heights. This indicates that the two
parameters introduced to model the dissipated power in practice take values
that lead to a flow rate compatible with the Hagen’s equation. It is worth
mentioning that the contribution due to the second term in Eq. (5), which
corresponds to the dissipation in the converging zone of the flow, is small
but not negligible (see dot-dashed line in Fig. 4).

###  VII.2 The bulk density and the prefactor $C$

As we mentioned in the introduction, Hagen’s prediction for the flow rate is
$Q=-\dot{M}=(\pi/4)\rho_{\textrm{o}}\sqrt{g}D_{\textrm{o}}^{5/2}$. Therefore,
by measuring the packing fraction $\rho_{\textrm{o}}$ at the opening and
$\rho_{\textrm{b}}$ in the bulk of the silo one can provide an estimate for
the constant $C$ since, according to Hagen’s analysis,
$\pi\rho_{\textrm{o}}/4=C\rho_{\textrm{b}}$. During discharge,
$\rho_{\textrm{o}}\approx\rho_{\textrm{b}}/2$ (see Appendix C for details),
implying that $C\approx 0.39$. Surprisingly, this estimate is more than $30\%$
below the value obtained by fitting experiments by various author
($C_{\textrm{experim}}\approx 0.58$) NeddermanBook .

In Eq. (V) the non-dimensional prefactor (i.e., $\pi\sqrt{2}/4$) differs from
the Hagen’s prefactor. If we use $\rho_{\textrm{o}}\approx\rho_{\textrm{b}}/2$
(see Appendix C) to replace the packing fraction at the orifice by the bulk
packing fraction in Eq. (V), the non-dimensional prefactor becomes
$\frac{\pi\sqrt{2}}{8}\approx 0.56$, which can be compared with the constant
$C$ in the Beverloo rule. This value deviates less than $4\%$ from the known
value $C_{\textrm{experim}}\approx 0.58$. Hence our Eq. (V) seems to be much
more accurate in predicting the prefactor that the Hagen equation.

It is important to note that the value of $C$ is determined by the expression
for the power lost through the orifice [see Eq. (34) and compare with Eqs. (7)
and (V)]. A key role is played by the area of the cross section of the orifice
$A_{\rm o}$. For a non-circular cross section, the factor $C$ will take a
different value. For example, for a square orifice of side $L$ ($A_{\rm
o}=L^{2}$), the prefactor becomes $C=\sqrt{2}/2$ (bear in mind that we take
$\rho_{\textrm{o}}\approx\rho_{\textrm{b}}/2$). This is consistent with the
experimental results from Ref. fowler1959 .

Based in the previous discussion, in the rest of the paper we will set
$C=\frac{\pi\sqrt{2}}{8}$ when we use the Beverloo expression for a circular
orifice.

###  VII.3 The $-kd$ boundary effect correction

When calculating our energy contributions in Appendix B, we assumed that the
system can be considered as a continuum while calculating
$\dot{W}_{\textrm{out}}$. We calculated the mass flow rate as
$-\dot{M}(t)=\rho_{\textrm{o}}A_{\textrm{o}}v_{\textrm{out}}$ without
considering that the orifice has a size only a few times the size of one
grain. This introduces boundary corrections that we did not take into account.

The simplistic correction, in the style of Beverloo, done by replacing the
orifice diameter $D_{\textrm{o}}$ by an effective smaller diameter
$D_{\textrm{o}}-kd$ has been questioned mankoc2007 ; janda2012 . However, this
is a simple way to incorporate this effect and we will use this in what
follows. For spherical grains, most authors conducting either experiments or
simulations in 3D indicate that $k\approx 1.4$. However, such fitting value
corresponds to the parameter $C$ fitted to $0.58$. As discussed in the
previous section, we now have a theoretical basis for setting
$C=\frac{\pi\sqrt{2}}{8}$. If we do this and fit the only remaining parameter
$k$ we find, for our DEM data (see Fig. 5) that $k=1.72\pm 0.02$. This is the
value we will use for the rest of the paper. Notice that with $C$ and $k$
selected in this manner the Beverloo equation and the asymptotic expression
(18) for the current theory coincide.

Figure 5: Mass flow rate as a function of $(D/d-k)^{5/2}$. Symbols correspond
to the DEM data. The lines correspond to both the Beverloo equation (yellow
solid) and the asymptotic expression (black dashed) for the present work (Eq.
(18) plus the boundary effect correction) setting $C=\frac{\pi\sqrt{2}}{8}$
and fitting $k=1.72\pm 0.02$. For these simulations $D_{\rm s}=30d$ and
$\mu=0.5$ for the particle–particle and the particle–wall interactions.

##  VIII Evolution of the flow rate during discharge

Equation (V) closed by Eqs. (14) and (VI.1) do not need to be solved
numerically as a function of time since a parametric plot of $-\dot{M}(t)$ as
a function of $M(t)$ can be directly obtained. In Fig. 6(a) we plot the flow
rate during discharge. We also include $M(t)$ without scaling in Fig. 6(b) as
a reference. As we can see, for most part of the discharge, Eq. (V) predicts a
constant flow rate, in agreement with our DEM simulations and experimental
observations. However, the model predicts an early drop in the flow rate well
before this is observed in the simulations. We believe this discrepancy is
connected with the poor prediction of the pressure contribution
$\langle\sigma_{rr}\rangle$ for the final stages of the discharge (see Fig.
3(a)). We expect that new developments on the estimation for the internal
pressure in silos will lead to an immediate improvement of the prediction of
Eq. (V) for the evolution of the flow rate.

Despite the shortenings of the prediction in the final stages of the
discharge, Eq. (V) provides, to our knowledge, the first mean to calculate the
evolution of the discharge. This opens new possibilities to study problems
such as forced discharges that show an non-constant flow rate madrid2018 .

Figure 6: (a) Mass flow rate as a function of the total mass $M$ in the silo
during the discharge. The flow rate is scaled by the value from the Beverloo
equation (i.e., the asymptotic flow rate for a tall column). The solid line
corresponds to the prediction by Eq. (V) using $\mu(I)=0.2$ and $\alpha=2.5$.
Symbols correspond to DEM data for different $D_{\rm o}$ and $D_{\rm s}$.
Unless otherwise stated $D_{\rm s}=30d$. The friction coefficient is set to
$\mu=0.5$ for the particle–particle and the particle–wall interactions. (b)
Mass in the silo as a function of time for the same simulations as in part
(a). Note that for the narrow silo we used a lower initial mass.

##  IX Predictions for low particle–particle friction

In the previous sections, we have observed an excellent agreement between the
Beverloo equation and our theory by using $\mu(I)=0.2$. This is possible
because $\mu(I)$ is not very sensitive to the material properties of the
grains if $\mu>0.4$ for the particle–particle interaction dacruz2005 . This is
why the a wide range of materials can be fitted to the Beverloo equation with
a single value for the constant $C$. However, if $\mu<0.4$ the effective
friction $\mu(I)$ starts to drop. Therefore, we must expect that the flow rate
will depend on the particle–particle friction if $\mu<0.4$. Some authors have
indeed reported that for low friction the flow rate is higher kondic2014 .
This effect is not accounted for in the Beverloo equation except for the fact
that the constant $C$ can be fitted to a new value. However, Eq. (V) can
predict the effect of lowering $\mu$ without changing neither the non-
dimensional prefactor nor $k$ (we recall that $k$ is introduced when
correcting $D_{\rm o}$ by $D_{\rm o}-kd$). The correction is obtained by
tunning $\mu(I)$, which has a clear physical interpretation from the
$\mu(I)$-rheology.

In Fig. 7(a), we plot the mass flow rate during discharge for different values
of $\mu$ used in our DEM simulations. We scaled the flow rate by the one
predicted by Beverloo in the previous sections, which holds valid for large
friction ($\mu>0.5$). The solid lines correspond to Eq. (V), where we have
used as a fitting parameter the value of $\mu(I)$. For $\mu>0.4$ we use
$\mu(I)=0.2$ as in the previous sections. For $\mu<0.4$ we set $\mu(I)$ to
lower values, while we keep $\alpha^{-1}=\tan\phi$ and $k=1.72$ as in the
previous section. The actual values of $\mu(I)$ used are shown in Fig. 7(b).
For reference, we also include $M(t)$ without scaling in Fig. 7(c).

The values of $\mu(I)$ are difficult to predict. These depend, for example, on
the geometry. In the quasistatic limit for plane shear, while $\mu(I)=0.28$ in
2D dacruz2005  it rises to $0.36$ in 3D azema2014 . It is difficult to provide
an independent estimate for $\mu(I)$ as a function of $\mu$ for the particular
case of a cylindrical silo. However, Eq. (V) provides with the new insight
that the mass flow rate scales as $\sqrt{\mu(I)}$.

Figure 7: (a) Mass flow rate as a function of mass in the silo for DEM
simulations for different friction coefficients $\mu$. We set $D_{\rm o}=10d$
and $D_{\rm s}=30d$. The particle–particle and particle–wall values of $\mu$
are made equal. The flow rate is scaled by the Beverloo equation. The solid
lines correspond to the prediction of Eq. (V) where $\mu(I)$ is set to $0.2$
if $\mu>0.4$ and to lower values if $\mu<0.4$. (b) Values of $\mu(I)$ as a
function of $\mu$ used in part (a). The solid line in part (b) is only a guide
to the eye. (c) Mass in the silo as a function of time for the same
simulations as in part (a).

##  X Conclusions

We have used the work–energy theorem to derive an expression for the mass flow
rate of a discharging silo as a function of the mass inside the silo [see Eq.
(V)]. For wide silos and stiff grains, we have shown that the changes in
kinetic energy and elastic energy can be neglected. We have used the concepts
of the $\mu(I)$-rheology to calculate the power dissipated during discharge.
The asymptotic limit of Eq. (V) resembles the well known Beverloo rule.
Interestingly, the non-dimensional prefactor that we predict is within $4\%$
of the experimental fitted values. Besides, we have shown that Eq. (V)
provides a mean to explain the higher flow rates observed in low-friction
materials.

It would be important to test the limitations of this approach when different
conditions are used such as: two-dimensional silos; silos, orifices and
particles with different shapes; hoppers; bumpy walls; use of an overweight;
use of soft and/or deformable grains; etc. In particular, recent simulations
with spherocylinders hidalgo2018  have shown that $\mu(I)$ depends on aspect
ratio, which should have an impact on the flow rate. It is also interesting
the potential application to suspensions and submerged grains passing through
constrictions since there are recent developments that indicate that the
$\mu(I)$-rheology is suitable to describe the flow in these systems boyer2011
; Houssais2016 . Forced silo discharges using overweights have shown non-
constant flow rates during discharge madrid2018 ; madrid2019 . An extension of
Eq. (V) to forced flows may be suitable to model the discharge under such
conditions.

Finally, we suggest a potential extension for this theory. The
$\mu(I)$-rheology is a local approximation that does not account for the
effects of particle size. The use of non-local approximations (see for example
kamrin2012 ) may help to obtain a correction to account for the “empty
annulus” effect and so avoid the simplistic approximation introduced to
consider the boundary effects at the orifice.

##  Appendix A Interactions for the DEM simulations

We used LIGGGHTS liggghts  to calculate the trajectories of each particle by
integrating the Newton-Euler equations poschel2005 . The equations of motion
for $N$ grains are solved via a velocity-Verlet algorithm by advancing in
small time intervals $\Delta t$. We used the Hertz particle–particle
interaction model with Coulomb criterion poschel2005 . In this model, the
normal component of the contact force between two grains $i$ and $j$ is
defined as

$F_{\rm n}=k_{\rm n}\delta_{\rm n}^{3/2}-\gamma_{\rm n}\delta_{\rm
n}^{1/4}\dot{\delta}_{\rm n},$ (19)

where

$k_{\rm n}=\frac{\sqrt{d}}{3}\frac{Y}{(1-\nu^{2})}$ (20)

is the elastic constant for normal contacts and $\delta_{\rm n}$ is the
overlap in the normal direction between particles $i$ and $j$. $Y$ is the
Young’s modulus, $\nu$ the Poisson ratio and $d$ the diameter of the
particles. The damping constant for normal contact is

$\gamma_{\rm
n}=\sqrt{\frac{5}{6}}\frac{\ln(e)}{\sqrt{\ln^{2}(e)+\pi^{2}}}\sqrt{\frac{Ym}{1-\nu^{2}}}\sqrt[4]{d},$
(21)

which results from the solution of the Herztian spring-dashpot model
antypov2011 . Here, $m$ is the grain mass and $e$ is de coefficient of
restitution wich is independent of the velocity tsuji1993 . $\dot{\delta}_{\rm
n}$ is the normal relative velocity of the particles $i$ and $j$.

The tangential component of the contact force between particles $i$ and $j$ is

$F_{\rm t}=-\textrm{sign}(v_{\rm t})\textrm{min}(|k_{\rm t}\delta_{\rm
t}\delta_{\rm n}^{1/2}-\gamma_{\rm t}\dot{\delta}_{\rm t}\delta_{\rm
n}^{1/4}|,\mu F_{\rm n}),$ (22)

being $v_{\rm t}$ the tangential relative velocity of the spheres at the point
of contact which takes into accounts the relative velocity of the centers of
the spheres and their rotation. The elastic constant for tangential contacts
is

$k_{\rm t}=\frac{Y\sqrt{d}}{(2-\nu)(1+\nu)}.$ (23)

The damping constant for tangential contact is

$\gamma_{\rm
t}=-2\sqrt{\frac{5}{6}}\frac{\ln(e)}{\sqrt{\ln^{2}(e)+\pi^{2}}}\sqrt{\frac{Ym}{2(2-\nu)(1+\nu)}}\sqrt[4]{d}.$
(24)

$F_{\rm t}$ is limited by Coulomb friction, being $\mu$ the friction
coefficient poschel2005 . In this model, the static and dynamic friction
coefficients coincide. The tangential displacement $\delta_{\rm t}$, which
depends on the history of the contact, is calculated as

$\delta_{\rm t}(t)=\int^{t}_{t_{\rm c}}v_{\rm t}(t^{\prime})dt^{\prime},$ (25)

where $t_{\rm c}$ is the time at which the contact begins.

To calculate the dissipated power in DEM simulations there exist two basic
approaches as described below.

_Approach 1._ At each time step the work done by the non-conservative terms in
the contact forces is calculated and this is saved in a cumulative variable.
This is done by taking the dot product between the displacement
$\mathbf{\delta r}$ of each contact and the non-conservative part of the
contact force in that time step [i.e., the force vector resulting from
combining the second term in Eq. (19) and either, the second term in the
absolute value of Eq. (22) or $\mu F_{\rm n}$, whichever applies]. This
approach has the advantage that one can accumulate separately different
contributions such as the tangential contribution [frictional, Eq. (22)] and
the normal contribution [inelastic collision, Eq. (19)]. However, this is CPU
demanding since the operation has to be done at every time step of the
simulation.

_Approach 2._ At two different arbitrary times (which may be separated by many
time steps) one calculates the potential energy of all conservatives forces
and the total kinetic energy (including rotations). The conservative forces
include gravity and the conservative terms of the contact forces [i.e., the
first term in Eq. (19) and the first term in the absolute value of Eq. (22)].
These potential energies can be calculated at any time since the conservative
forces depend only on the current positions of the particles (not the actual
trajectories nor velocities). Therefore, one does not need to make this
calculation at every time step to track trajectories. The difference of total
energy (potential + kinetic) between the two times under consideration
corresponds to the energy dissipated by the non-conservative terms in the
contact interactions. Unfortunately, this approach will not provide detailed
information on the contribution of each dissipation mode (friction and normal
collisions).

We have used _Approach 2_ to obtain the dissipated energy. This prevents us
from having access to the detail of how much energy is lost by friction and
how much by normal collisions. However, this is a very efficient method and
suffices for the purposes of the current study.

##  Appendix B Contributions to the work–energy theorem

In Fig. 8, we plot the different contributions to the power injected or
extracted from the silo obtained via a DEM simulation of the silo discharge.
We have run simulations for various silos diameters, particle–particle
interaction parameters and orifice sizes. All cases studied display the same
trends as the ones shown in this sample simulation.

###  B.1 Internal kinetic energy ($\dot{K}_{\textrm{in}}$)

As we can see in Fig. 8, the contribution of $\dot{K}_{\textrm{in}}$ is one
order of magnitude smaller than $\dot{W}_{\rm out}$ and three orders of
magnitude smaller than $\dot{W}_{\rm g}$. Therefore, we can neglect
$\dot{K}_{\textrm{in}}$. To provide a partial explanation for this observation
consider the kinetic energy of the grains inside the silo madrid2017b

$K_{\textrm{in}}(t)=\frac{1}{2}\sum_{i=1}m\bm{v}_{i}^{2}(t).$ (26)

The sum runs over all particles inside the silo at time $t$ and
$\bm{v}_{i}(t)$ is the velocity of particle $i$. This can be expressed in
terms of the center of mass velocity $v_{\textrm{cm}}$ and the “granular
temperature” as
$K_{\textrm{in}}(t)=K_{\textrm{in}}^{\textrm{cm}}(t)+K_{\textrm{in}}^{\textrm{temp}}(t)$,
being

$K_{\textrm{in}}^{\textrm{cm}}(t)=\frac{1}{2}M(t)v_{\textrm{cm}}^{2}(t),$ (27)

and

$K_{\textrm{in}}^{\textrm{temp}}(t)=\frac{1}{2}\sum_{i=1}^{N(t)}m_{i}[\bm{v}_{i}(t)-\bm{v}_{\textrm{cm}}(t)]^{2}.$
(28)

We will disregard the “temperature” term and focus on the center of mass. We
assume here that the $x$ and $y$ components of $\bm{v}_{\textrm{cm}}$ are null
and therefore $|\bm{v}_{\textrm{cm}}|=v_{\textrm{cm}}=-\dot{z}_{\textrm{cm}}$.
Hence, the rate of change of the kinetic energy is

Figure 8: Contributions to the work–energy balance during the discharge of a
silo as a function of time (DEM simulations). Power injected by gravity
$\dot{W}_{\rm g}$, dissipated power $\dot{W}_{\rm d}$, power loss through the
orifice $\dot{W}_{\rm out}$, elastic power $\dot{W}_{\rm el}$ and rate of
change of the internal kinetic energy $\dot{K}_{\rm in}$. Data corresponds to
a silo with $D_{\rm s}=30d$ and $D_{\rm o}=6d$ while the particle–particle and
particle–wall friction coefficients are set to $\mu=0.5$.

$\displaystyle\dot{K}_{\textrm{in}}(t)$
$\displaystyle\approx\dot{K}_{\textrm{in}}^{\textrm{cm}}(t)$
$\displaystyle=M(t)v_{\textrm{cm}}(t)\dot{v}_{\textrm{cm}}(t)+\frac{1}{2}v_{\textrm{cm}}^{2}(t)\dot{M}(t),$
(29)

which can be written, using Eqs. (2) and (3), as

$\dot{K}_{\textrm{in}}(t)\approx\frac{1}{4\rho^{2}_{\textrm{b}}A^{2}_{\textrm{s}}}\left[M(t)\dot{M}(t)\ddot{M}(t)+\frac{1}{2}\dot{M}^{3}(t)\right].$
(30)

As we can see, $\dot{K}_{\textrm{in}}$ decays as $A_{\textrm{s}}^{-2}$.
Therefore, $\dot{K}_{\textrm{in}}$ will be small for wide silos. We will see
below that other terms in the energy balance decrease as
$A_{\textrm{s}}^{-1}$.

###  B.2 Gravitational energy ($\dot{W}_{\textrm{g}}$)

The gravitational potential energy of the particles inside the silo is

$U_{\textrm{g}}(t)=M(t)gz_{\textrm{cm}}(t)=\frac{gM^{2}(t)}{2\rho_{\textrm{b}}A_{\textrm{s}}}.$
(31)

Where we have used Eq. (2). Therefore, the power injected by the action of
gravity is

$\dot{W}_{\textrm{g}}(t)=-\dot{U}_{\textrm{g}}(t)=-\frac{g}{\rho_{\textrm{b}}A_{\textrm{s}}}M(t)\dot{M}(t).$
(32)

As we can see, $\dot{W}_{\textrm{g}}$ scales with $A_{\textrm{s}}^{-1}$ in
contrast to the faster decay displayed by $\dot{K}_{\textrm{in}}$. We recall
here that $\dot{M}(t)$ is negative; therefore $\dot{W}_{\textrm{g}}(t)$ is
positive. It is worth mentioning that $\dot{W}_{\textrm{g}}$ is mostly
dissipated and goes into $\dot{W}_{\rm d}$ (see Fig. 8). A very small portion
of $\dot{W}_{\textrm{g}}$ goes into the draining grains
($\dot{W}_{\textrm{out}}$), as we discuss in the following subsection.

###  B.3 Discharge energy loss ($\dot{W}_{\textrm{out}}$)

While the system discharges, the particles that leave the system take away
some energy since their own kinetic energy is no longer part of the internal
energy of the silo. If during a time interval $dt$ at time $t$ the system
discharges a mass $dM=-\dot{M}(t)dt$ at velocity $v_{\rm out}$, then the
kinetic energy removed per unit time is

$\dot{K}_{\textrm{out}}(t)=-\frac{1}{2}v_{\textrm{out}}^{2}\dot{M}(t).$ (33)

By definition, $-\dot{M}(t)=\rho_{\textrm{o}}A_{\textrm{o}}v_{\textrm{out}}$.
Therefore,

$\dot{W}_{\textrm{out}}=\dot{K}_{\textrm{out}}(t)=-\frac{\dot{M}^{3}(t)}{2\rho^{2}_{\textrm{o}}A^{2}_{\textrm{o}}}.$
(34)

The power lost by discharge does not depend on the silo cross section but on
the orifice cross section $A_{\textrm{o}}$. Then again,
$\dot{W}_{\textrm{out}}$ is positive since $\dot{M}(t)$ is negative.

It is important to note that the definition
$-\dot{M}(t)=\rho_{\textrm{o}}A_{\textrm{o}}v_{\textrm{out}}$ is not strictly
correct since the velocity and the density at the orifice are not homogeneous
(see Appendix C for details) janda2012 .

###  B.4 Elastic energy ($\dot{W}_{\textrm{el}}$)

If the grains are stiff, the variation in the elastic energy at the contacts
are expected to be small. $\dot{W}_{\textrm{el}}$ corresponds to the rate of
change of the conservative component of the contact forces
$\dot{E}_{\textrm{el}}$. In Fig. 8, we show the rate of change in elastic
energy during a silo discharge in a DEM simulation. As we can see, this term
is of the order of the change in kinetic energy $\dot{K}_{\rm in}$ and can be
neglected in comparison with $\dot{W}_{\rm g}$ and $\dot{W}_{\rm out}$. Of
course, this may be inadequate for very soft grains. For our purposes, we
simply disregard this contribution in the analysis. However, this can be
eventually included by using the stress based expression for the elastic
energy density in the limit of small deformations landau

$E_{\textrm{el}}=A_{\textbf{s}}z_{\text{cm}}\bm{\sigma}\bm{\epsilon},$ (35)

being $\bm{\sigma}$ the stress tensor and $\bm{\epsilon}$ the strain tensor.

##  Appendix C Estimation of $\rho_{\rm o}$

Since the packing fraction and the particle velocities are not homogeneous
across the orifice, the correct flow rate has to be calculated as janda2012

$Q=\int_{0}^{D_{\rm o}/2}{v_{\rm o}(r)\rho_{\rm o}(r)2\pi rdr},$ (36)

where $v_{\rm o}(r)$, and $\rho_{\rm o}(r)$ are the vertical velocity and
density profiles at the horizontal plane of the orifice as a function of the
distance $r$ to the center of the orifice.

If the profiles where flat, one can write

$Q=\rho_{\rm o}v_{\rm o}A_{\rm o},$ (37)

where $A_{\rm o}$ is the cross section of the orifice and $\rho_{\rm o}$ and
$v_{\rm o}$ are the mean density and mean vertical velocity at the orifice,
i.e.,

$\displaystyle\rho_{\rm o}=A_{\rm o}^{-1}\int_{0}^{D_{\rm o}/2}{\rho_{\rm
o}(r)2\pi rdr},$ $\displaystyle v_{\rm o}=A_{\rm o}^{-1}\int_{0}^{D_{\rm
o}/2}{v_{\rm o}(r)2\pi rdr}.$ (38)

This is in general invalid (for non-flat profiles) because the integral of the
product $v_{\rm o}(r)\rho_{\rm o}(r)$ is different from the product of the
integrals of $v_{\rm o}(r)$ and $\rho_{\rm o}(r)$. It is worth mentioning that
if $\rho_{\rm o}(r)$ does not depend on $r$, then Eq. (37) is valid if $v_{\rm
o}$ is defined as in Eq. (C). This is why expression (37) is valid for
incompressible fluids. However, for granular materials, the density across the
orifice is not constant and Eq. (37) is actually a poor approximation if
$\rho_{\rm o}$ and $v_{\rm o}$ are defined as in (C).

Therefore, the definitions for $v_{\rm o}$ and $\rho_{\rm o}$ in the
simplified Eq. (37) need to be replaced by “effective values” rather than
“mean values”. As an example, we can set $v_{\rm o}$ to the mean value of
$v_{\rm o}(r)$ but then chose $\rho_{\rm o}$ to yield the correct flow rate.
We have done this for experiments and simulations from the literature and
found that setting $v_{\rm o}$ to the mean value of $v_{\rm o}(r)$ leads to an
effective value of $\rho_{\rm o}\sim\rho_{\rm b}/2$. The mean value for the
packing fraction at the orifice is in fact lower (about $\rho_{\rm b}/3$). As
an example, in our simulations for $D_{\rm o}=6d$, the mass flow rate is
$Q=0.00312$ kg/s. The corresponding vertical velocity at the orifice (averaged
over a cylindrical region that fits exactly in the orifice and has a height of
$d$) is $v_{\rm o}=0.1353$ m/s. Therefore, from Eq. (37) above, we find that
the effective apparent density at the orifice is $\rho_{\rm o}=815$ kg/m3.
Since our bulk $\rho_{\rm b}=1475$ kg/m3, then $\rho_{\rm o}/\rho_{\rm
b}=0.55$.

###### Acknowledgements.

This work has been supported by ANPCyT (Argentina) through grant
PICT-2012-2155, Universidad Tecnológica Nacional (Argentina) through grant
PID-MA0FALP0002184 and PID IFI1871, Centro Argentino Francés de Ciencias de la
Ingeniería (CAFCI, Argentina-Francia) and FONACIT through grant 2015000072
(INVUNI2013-1563) Universidad Simón Bolívar (Venezuela).

## References

  * (1) G. H. L. Hagen, Bericht uber die zur Bekanntmachung geeigneten Verhandlungen der K oniglich Preussischen Akademie der Wissenschaften zu Berlin. 35–42 (1852).
  * (2) B. P. Tighe and M. Sperl, Gran. Matt. 9, 141 (2007)
  * (3) R. L. Brown and J. C. Richards, Principles of Powder Mechanics, Pergamon Press, Oxford, (1970).
  * (4) W. Beverloo, H. Leniger and J. Van de Velde, Chem. Eng. Sci. 15, 260 (1961).
  * (5) R. M. Nedderman, Statics and kinematics of granular materials (Cambridge University Press, Cambridge, 2005).
  * (6) M. A. Madrid, J. R. Darias and L. A. Pugnaloni, Forced flow of granular media: Breakdown of the Beverloo scaling, Europhys. Lett. 123, 14004 (2018).
  * (7) R. M. Nedderman, U. Tüzün, S.B. Savage and G.T. Houlsby, J. Chem. Eng. Sci. 37, 1597-1609 (1982)
  * (8) U. Tüzün, G. T. Houlsby, R. M. Nedderman and S. B. Savage, J. Chem. Eng. Sci.37, 1691-1709 (1982)
  * (9) S. B. Savage, R. M. Nedderman, U. Tüzün and G.T. Houlsby, J. Chem. Eng. Sci. 38, 189-195 (1983).
  * (10) G. H. Ristow, Physica A 235, 319 (1997).
  * (11) C Mankoc, A Janda, R Arévalo, J Pastor, I Zuriguel, A Garcimartín and D Maza, Granular Matter 9, 407 (2007).
  * (12) L Kondic, Simulations of two dimensional hopper flow, Granular Matter 16, 235 (2014).
  * (13) T J Wilson, et al., Granular discharge rate for submerged hoppers, Papers in Physics 6, 060009 (2014).
  * (14) E. Goldberg, C. M. Carlevaro and L. A. Pugnaloni, Papers in Physics 7, 070016 (2015).
  * (15) M. C. Villagrán Olivares, J. G. Benito, R. O. Uñac and A. M. Vidales, Towards a one parameter equation for a silo discharging model with inclined outlets, Powder Technology 336 265 (2018).
  * (16) S. M. Rubio-Largo, A. Janda, D. Maza, I. Zuriguel and R. C. Hidalgo, Phys. Rev. Lett. 114, 238002 (2015).
  * (17) I. R. McDougall and A. C. Evans, Rheologica Acta 4, 218 (1965).
  * (18) L. Staron, P.-Y- Lagrée S. Popinet, Phys. Fluids 24, 103301 (2012).
  * (19) F. da Cruz, S. Emam, M. Prochnow, J. N. Roux and F. Chevoir, Phys. Rev. E 72, 021309 (2005).
  * (20) C. Kloss, C. Goniva, A. Hager, S. Amberger and S. Pirker, Prog. Comput. Fluid Dynamics, 12, 140 (2012).
  * (21) E. Azéma and F. Radjai, Phys. Rev. Lett. 112, 078001 (2014).
  * (22) Note that we have selected the negative root since this is the physically meaningful root corresponding to negative $\dot{M}$.
  * (23) J. K. Walters, A theoretical analysis of stresses in silos with vertical walls, Chem. Eng. Sci. 28, 13 (1973).
  * (24) D. M. Walker, An approximate theory for pressures and arching in hoppers, Chem. Eng. Sci. 21, 975 (1966).
  * (25) R. T. Fowler and J. R. Glastonrury, Chem. Eng. Sci. 10, 150 (1959).
  * (26) A. Janda, I. Zuriguel and D. Maza, Phys. Rev. Lett. 108, 248001 (2012).
  * (27) R. C. Hidalgo, B. Szabó, K. Gillemot, T. Börzsönyi, T. Weinhart, Phys. Rev. F 3, 074301 (2018).
  * (28) F. Boyer, E. Guazzelli and O. Pouliquen, Phys. Rev. Lett. 107, 188301 (2011).
  * (29) M. Houssais, C. P. Ortiz, D. Durian and J. Jerolmack, Phys. Rev. E 94, 062609 (2016).
  * (30) M. A. Madrid and L. A. Pugnaloni, Granular Matt. 21, 76 (2019).
  * (31) K. Kamrin and G. Koval, Nonlocal constitutive relation for steady granular flow, Phys. Rev. Lett. 108, 178301 (2012).
  * (32) T. Pöschel and T. Schwager, Computational Granular Dynamics. Models and Algorithms, Springer, Berlin (2005).
  * (33) D. Antypov and J. Elliott, Europhys. Lett. 94, 50004 (2011).
  * (34) Y. Tsuji and T. Kawaguchi, T. Tanaka, Powder Technol. 77, 79 (1993).
  * (35) M. A. Madrid, J. R. Darias and L. A. Pugnaloni, EPJ Web of Conf. 140, 03041 (2017).
  * (36) L. D. Landau and E. M. Lifshitz, Theory of Elasticity (3rd ed.). Butterworth Heinemann, Oxford (1986).