1. Introduction
Understanding the biomechanical features and electromechanical response of biological tissues and biomedical components is of significance in many medical science and technology fields including implantable devices, tissue-engineered constructs, drug delivery systems and disease detection [
1,
2]. To increase the efficiency and durability of biomedical implantable devices such as knee joint replacement tools, cardioverter defibrillators, and dental implants, mechanical features including elasticity properties and vibrational characteristics have vital roles since they are highly associated with the overall performance of the device [
3,
4]. Biomechanical properties are also important in the design and manufacture of biomaterials used as the building blocks of artificial organs and tissue-engineering scaffolds [
5,
6]. Optimised mechanical properties lead to the smooth and robust integration of biomaterials into target tissues, reducing inflammation and enhancing the efficacy of the artificial organ or scaffold [
5,
7]. In drug delivery systems, acoustic analysis and fluid flow studies help researchers and clinicians gain a better understanding of therapeutic particle motion and improve drug delivery rates [
8,
9]. Moreover, biomechanical features are used as indicators of disease for detection purposes and in treatment monitoring to assess therapy effectiveness [
10,
11,
12].
In clinical applications, ultrasound transducers are commonly used to create waves within tissue and the patient’s body [
13]. Depending on the direction in which the ultrasound transducer is applied by the professional nurse or general practitioner, it can create in-plane, shear, and flexural waves [
14,
15]. If the transducer as the initial source of the waves is applied perpendicular to the human tissue, it leads to flexural waves; meanwhile, when the ultrasound transducer is utilised along the side of the tissue or human body, the in-plane wave propagation assumption is valid [
14,
15]. Each direction of wave propagation provides useful information about potential diseases or abnormalities. A combination of them is used in some applications for better spatial visualisation [
13,
15].
Elastography is a promising non-invasive biomedical imaging technique in which the biomechanical features of the target biological tissue are visualised on an image called an elastogram for disease detection or treatment efficacy assessments [
15,
16]. Depending on the source of the mechanical stimulation, there are several types of elastography imaging including ultrasound [
17], magnetic resonance [
18], and optical [
12]. In all of these applications, mechanical deformation is created using a safe and non-invasive approach and then the corresponding displacements per unit length (strain) are measured. The obtained stress and strain are utilised in a mathematical model to estimate the biomechanical features of the target tissue [
19]. Elastography imaging has been used in the rapid and quantitative detection of many diseases and abnormalities, especially those in which fibrosis, stiffness hardening, or/and softening are involved. These biomedical imaging applications include—but are not limited to—breast cancer [
20], liver fibrosis [
21], polycystic ovary syndrome [
22], and prostate cancer [
23].
Accurate and computationally fast mathematical models are essential in elastography imaging in order to calculate and plot intrinsic biomechanical features such as stiffness [
24], Poisson’s ratio [
25], and poroelastic [
26] and viscoelastic properties [
27]. To provide elastography mathematical models of deformation for clinical applications, different continuum mechanics-based models have been developed in recent years. Samani et al. [
28] introduced a finite element approach for studying the elasticity behaviour of human breast samples for cancer detection. They conducted finite element simulations and experimentally obtained the elasticity modulus of a wide list of different pathological breast tissues including both malignant and benign conditions. It was found that the presence of cancer in human breast tissue leads to a significant increase in the elasticity modulus, which could be clinically beneficial in the detection of breast cancer. Kheirkhah et al. [
29] used mathematical constraints obtained based on the principles of tissue mechanics for the regularization of displacements and strains in ultrasound elastography, incorporating factors such as tissue incompressibility and deformation compatibility. Poul et al. [
30] studied different rheological models of viscoelasticity for elastography imaging applications including the Kelvin–Voigt, standard linear solid, and Kelvin–Voigt fractional derivative models. They concluded that the Kelvin–Voigt fractional derivative model performs consistently in terms of stress relaxation. Moreover, hyperelastic constitutive models [
31] and porohyperviscoelastic approaches [
32] have been developed for shear-wave elastography and have been investigated for tissue-mimicking materials and human liver, respectively. A doublet mechanics model [
33], a convolutional neural network technique [
34], and an analytical poroelastic approach [
35] have also been proposed to improve tissue mechanical property visualisation. For more details about the state-of-the-art methods and techniques used to obtain the biomechanical features of biological tissues, an interested reader is referred to a recent review paper by Alekya et al. [
36].
In recent years, several higher-order biomechanical models have been developed to study the deformation behaviour of biological components and tissues under mechanical forces. Xiang and Liew [
37] introduced a highly accurate computational platform to analyse the transverse compression of protein microtubules using a higher-order Cauchy–Born law. Moreover, Xiang et al. [
38] developed a Cauchy–Born rule of third order in conjunction with an element-free computational framework to simulate the biomechanical behaviour of microtubules. More recently, Dwairy et al. [
39] applied the biphasic theory to predict interstitial fluid pressure and stress within solid tumours. They found that mechanical stress in solid tumours highly depends on the nonlinearity and constitutive model of deformation while the interstitial fluid pressure is relatively independent of the fundamental constitutive equation.
However, the current mathematical models that describe the mechanical properties of biological tissues do not account for the fundamental differences in tissue softening and hardening observed at different scales. Recent experimental research studies have revealed that the biomechanical features of biological tissues depend on the scale (i.e., size) at which the investigation is conducted [
40]. Fuhs et al. studied the stiffness of breast cancer cells as well as the elasticity modulus of human breast tissues and found that breast cancer cells are softer than healthy cells while cancer tissue is more rigid compared to the healthy tissue [
40]. This paradox cannot be explained in the context of classical continuum mechanics models as they are limited to large-scale dimensions and can only predict stiffness hardening at the tissue level [
41].
To the best of our knowledge, no scale-dependent mathematical models have been proposed to solve this paradox (i.e., the rigid-tumours–soft-cells paradox). To address this, a scale-dependent modified continuum model incorporating stress nonlocality and strain gradient effects for the in-plane wave propagation analysis of biological tissues is developed. Particular attention is paid to human breast lesions to visualise the results of the scale-dependent model. The higher-order nonlocal strain gradient continuum mechanics theory, in conjunction with Hamilton’s law, is employed to derive the differential equations of the breast tissue. An analytical solution is obtained for the in-plane wave propagation in biological tissues. A deep learning model made of one flatten layer, several hidden layers, and a few dense layers is also developed and trained to improve the accuracy and allow for extra flexibility to better fit the clinical data. In addition, other machine learning models such as linear, Ridge, Lasso, Random Forest, and ElasticNet regression models are developed for comparison purposes. To verify the accuracy of the analytical approach, a finite element method is also implemented to extract the in-plane wave characteristics of healthy and diseased breast tissue including adipose and fibroglandular tissue, fibroadenoma, invasive ductal carcinoma, infiltrating lobular carcinoma, ductal carcinoma in situ, fibrocystic disease, fat necrosis, and invasive mucinous carcinoma.
Figure 1 shows the flowchart of this research study with the four main steps, namely, (1) problem specification, (2) development of the higher-order nonlocal model, (3) machine learning modelling, and (4) finite element verification study. The new higher-order mathematical formulation integrated with machine learning allows for the introduction of three new parameters in biological tissue modelling: (1) zero-order nonlocal parameter, (2) first-order nonlocal parameter, and (3) strain gradient parameter. The combination of these parameters enables the description of both stiffness softening and hardening, which have been observed at cell and tissue levels, respectively, and solves the rigid-tumours–soft-cells paradox. The proposed mathematical model and simulation data have applications in biomedical elastography imaging for the early and accurate detection of breast cancer.
2. Theoretical Modelling
2.1. A Higher-Order Nonlocal Model of Biological Tissue
The main novel aspect of this research work is the development of the first scale-dependent, higher-order nonlocal model for analysing in-plane elastic waves propagated within human breast tissue. In addition, the higher-order model of nonlocality is integrated with a deep learning neural network model, which is later described in
Section 2.3.
To the best of our knowledge, this is the first refined combination of a higher-order nonlocal strain gradient model and a deep learning model for the extraction of the wave propagation characteristics of biological samples at different scales.
In Eringen’s nonlocal elasticity theory [
42], stress at a given point of the tissue is a function of the strain at all points of the domain (stress nonlocality). This is the basic assumption of the nonlocal continuum mechanics [
42,
43] and is initially stated in the form of an integro-partial differential equation between stress and strain. However, as this integro-partial differential equation is challenging and difficult to use for modelling purposes, Eringen introduced a simpler differential form and proved that it could incorporate stress nonlocality in the mechanical analysis [
42,
43]. Moreover, he found that experimental data and lattice dynamics simulations supported this differential form very nicely [
42,
43]. In this section, a more advanced version of the nonlocal continuum mechanics known as the higher-order nonlocal strain gradient theory [
44] is used, which can incorporate not only the zeroth-order stress nonlocality but also the first-order one.
The classical elasticity theory assumes that the stress at a given point is only a function of the strain at that point. This does not allow for the incorporation of stress nonlocality, which has been widely reported to be associated with scale effects and is responsible for stiffness softening. Moreover, in the classical local models, the influence of the strain gradient is neglected, which is a complementary scale parameter associated with stiffness hardening at small scales. The combination of these effects within the framework of a newly introduced mathematical model solves the rigid-tumour–soft-cell paradox as it allows for both a reduction and increase in the stiffness due to the cancer’s formation and growth. Based on the higher-order version of the nonlocal strain gradient model (HNSM), the stress is related to strain as [
44,
45,
46]
where
and
where
,
,
, and
are the stress tensor, elastic property matrix, strain, and strain gradient parameter, respectively. Furthermore, the coefficients in Equation (2) are given by
and
, which are, respectively, the first and second nonlocal stress parameters [
46,
47,
48,
49].
,
, and
are the combined, first-order, and zeroth-order nonlocal operators, respectively. When the zeroth-order nonlocal operator is equal to the first-order one (i.e.,
), the HNSM is reduced to the lower-order nonlocal strain gradient model (LNSM) as in [
50,
51].
On the other hand, the classical nonlocal model (CLNM) is recovered when the nonlocal operators are identical and the strain gradient operator is zero, namely,
The CLNM has been successfully utilised to analyse the biomechanics of protein microtubules [
52,
53], providing evidence on the scale-dependency behaviour of biological components and the importance of nonlocal effects [
54,
55]. If only the strain gradient effects are included, both nonlocal operators are neglected, and the strain gradient operator is maintained (
), leading to the following constitutive equation of biological tissue.
Equation (6) describes the fundamental stress–strain relation of the classical strain gradient model (CLSM) [
56,
57,
58]. Moreover, the constitutive equation of the HNSM can be reduced to the classical elasticity model (CLEM) when
which is mathematically described as:
In this section, a higher-order model of elasticity with strain gradient and stress nonlocality effects is developed for the in-plane wave propagations in rectangular-shaped, thin biological tissue sections (
Figure 2). Using Equations (1) and (2) in the two-dimensional Cartesian coordinate system, the stress components are linked to the strain components as
and
are the Poisson’s ratio and elasticity modulus of the biological tissue.
For small deformation induced by the in-plane wave propagation, the strain components are given by
where
u and
v are the in-plane displacements along the
x and
y axes, respectively. The force stress resultants can be obtained from Equations (8) and (9) as
where the force stress resultants (
Tij) are defined by
Here,
h denotes the thickness of the tissue section.
Aij are the in-plane stiffness components of the tissue, which are calculated using
where
G12 is the shear elasticity modulus of the biological tissue. Using Hamilton’s law for motions in the
x-
y plane, one can obtain
where
Here,
is the biological tissue mass density. It is worth mentioning that mass density is entirely different from mammographic breast density. In the context of breast tissue, mammographic density refers to the relative amount of fibroglandular tissue to fat tissue, which is an independent risk factor for breast cancer and is also a significant problem in mammography imaging as it masks cancer [
59,
60,
61]. By contrast, mass density in Equation (14) refers to the mass per unit volume of the tissue, which is measured in terms of kg/m
3. Using Equations (10) and (13), the governing differential equations for predicting the deformation behaviour of biological tissues under in-plane wave propagation are derived, which are given in
Appendix A as Equations (A1) and (A2). It should be noted that, in the present analysis, the effect of energy loss is not considered. In terms of the mechanical response of breast tissues, it has been shown that, after an initial transient period, the system reaches a steady state, and the assumption of negligible energy loss becomes valid [
35].
To determine the in-plane wave propagation characteristics of the breast tissue, the solution of Equations (A1) and (A2) can be assumed as
where
V and
U represent the in-plane vibration amplitudes along the y and x axes, respectively. The in-plane wave numbers along the
y and
x axes are
ky and
kx, respectively. Moreover,
is the frequency of vibrations induced by in-plane waves. Substituting Equation (15) into the in-plane motion equations, namely, Equations (A1) and (A2), leads to
where
,
, and
represent the total stiffness matrix, total mass matrix, and in-plane vibration frequency, respectively. The square root of the eigenvalues of Equation (16) are the frequencies of the in-plane vibrations within the tissue induced by in-plane waves. The elements of matrix
and
are listed in
Appendix A.
2.2. Finite Element Simulation
To verify the accuracy of the higher-order nonlocal modelling, finite element simulations have been conducted using COMSOL Multiphysics Simulation software version 5.5. The Solid Mechanics (solid) module in conjunction with a linear elastic material model was used for the analysis. The general structural analysis of three-dimensional bodies was adopted for a more accurate analysis. The eigenfrequency solver technique was the multifrontal massively parallel sparse direct (MUMPS), which is a powerful rapid numerical method for solving large linear systems using a parallel computational approach. Three different breast tissue samples were considered for the finite element study. First, fibroglandular tissue was tested as the main source of mammographic breast density, which causes a masking effect problem in mammography imaging [
62,
63]. The second case study has been performed on ductal carcinoma in situ (DCIS) samples, which are widely known as the early form of human breast cancer [
64]. DCIS cancers are noninvasive, meaning that they have not invaded surrounding breast tissue or spread into other parts of the body. Detecting DCIS enables more effective and faster treatment, improving the survival rates of breast cancer patients. The third samples, which have been considered for the finite element study, were fat necrosis as a benign condition that might be misdiagnosed as cancerous tumours. They are commonly formed after a trauma or injury [
65]. The finite element simulations were conducted for the first ten mode shapes of tissue samples. A rectangular-shaped geometry with a thickness of 1 cm, length of 10 cm, and width of 10 cm was created in a three-dimensional coordinate system as the tissue section. Scale effects including the stress nonlocality and strain gradients were neglected as the available commercial finite element software packages lack the ability to simulate these effects. After the geometry creation, the elastic properties such as Poisson’s ratio and Young’s modulus were assigned to the material (tissue sections) based on the available clinical data in the literature [
28]. The displacement and velocity components were initialised to zero at the start of the analysis, indicating a zero or at-rest initial condition. The boundary conditions of the samples were assumed to be the following: the displacement components along the
x and
z axes were zero at
x = 0 and
x = a, the displacement components in the
y and
z directions were set to zero at
y = 0 and
y = b, and a free type of boundary condition was adopted for the top and bottom surfaces of the rectangular sample. A combination of tetrahedra and triangles (element type) was used to create a smooth and well-distributed mesh pattern with an average extra-fine mesh size. An eigenfrequency technique was utilised as the solver, which is suitable for eigenvalue problems with eigenmodes.
2.3. Artificial Neural Network Simulation
To capture the different mechanical behaviours of breast lesions at different scales and include both stiffness hardening and softening, a machine learning model that can account for nonlocality, complexity, and nonlinear patterns is required. In addition, the model needs to have the capability of dealing with a limited number of samples and small datasets while preventing overfitting. Deep neural network models, especially with rectified linear unit (ReLU) activation functions, provide the ability to incorporate the complexity of the in-plane wave characteristics of breast lesions and their nonlocal deformation behaviour at small scales. Furthermore, neural network models allow for flexibility in the architecture by offering different hyperparameter adjustment options including the change in the number of hidden layers and neurons as well as the choice of activation function, optimiser, and learning rate. This flexibility in the fine-tuning process is essential for optimising the performance of the machine learning model and for preventing overfitting when the training dataset is relatively small, as it is in many clinical diagnosis applications. Overall, the benefits of neural network models include—but are not limited to—a fast training process, flexibility in fine-tuning, the ability to predict complexity, and the capability of mitigating overfitting.
TensorFlow, NumPy, and pandas were used to develop the deep neural network model, to work with multi-dimensional matrixes, and to create the database, respectively [
66,
67]. Keras, as a high-level application programming interface for conveniently working with TensorFlow 2.14.0 machine learning software, was adopted [
68]. Moreover, to divide the dataset into 70% training and 30% testing subsets, the train_test_split method of the model_selection module of the sklearn python library was used. The dataset used for training and testing contained 100 human samples diagnosed with DCIS. To build the deep learning model, several architectures were considered using the sequential method of Keras. Overall, the neural network was composed of one flatten layer, a couple of hidden dense layers, and an output layer as shown in
Figure 3a. In clinical applications, it is often challenging to find patients of all different pathological conditions, leading to relatively small datasets. Thus, a simple neural network can be adopted to prevent overfitting and ensure the generalisation capability of the model despite the restricted human breast samples. In addition, simpler neural network algorithms are computationally inexpensive and faster to implement as they require less training time as well as a faster hyperparameter adjustment process. A sequential deep neural network model of feedforward type was utilised in this study for prediction, which was made of several hidden layers. A rectified linear unit (ReLU) was employed as the activation function of dense layers, which allows for the incorporation of nonlinearity and complexity in real-world problems. An Adam optimiser with a learning rate of 0.001 and mean-square-error loss function was used to compile the machine learning model. Data scaling was performed prior to the training process to ensure consistency with the scale of all features. The training process was conducted for 500 epochs. During each epoch, the entire training dataset was passed forward and backwards within the artificial neural network. The rationale for choosing a relatively high epoch number is the limited sample size, which is usually the case in clinical studies where providing human tissue samples is not straightforward. This ensures that the neural network model has been sufficiently exposed to the available data for training and generalisation purposes.
Figure 3b shows a flowchart to visualise the different steps and processes involved in the machine learning modelling of the in-plane wave propagation analysis of human breast lesions. The model evaluation, which is essential for hyperparameter adjustment, was conducted using the mean square error of both training and test datasets. The learning rate, number of neurons within a dense layer, number of dense layers, activation function, number of epochs, and the optimiser technique have been adjusted during the hyperparameter adjustment process to obtain a highly accurate neural network model.
In this research work, three different architectures of neural networks were considered: (1) FL-3HDL(512,128,10): one flatten layer (FL), three hidden layers (HDL) with 512, 128, and 10 neurons as well as one dense layer with one neuron as output layer (5 layers and 651 neurons in total); (2) FL-2HDL(128,10): one FL, two HDLs with 512 and 10 neurons, as well as one dense layer with one neuron; (3) FL-1HDL(10): one FL, one HDL with 10 neurons, and one output dense layer. In all three different architectures, the activation function of the hidden and dense layers was set to the rectified linear unit (ReLU) to incorporate nonlocality and nonlinearity and to capture complex patterns within the data. The Adam optimiser with a learning rate of 0.001 was used to compile the machine learning model. The results of these three neural network models will be compared and discussed in the next section.
3. Results
To verify the analytical higher-order nonlocal strain gradient modelling, a finite element study was conducted on three different types of breast tissue, including healthy fibroglandular tissue (FGT), ductal carcinoma in situ (DCIS), and fat necrosis, which is a benign condition. The results (dimensionless phase velocities) are listed in
Table 1 for the first three mode numbers (
n = 1, 2, 3). The elasticity moduli of different human breast tissue were adopted from those reported by Samani et al. [
28]. They conducted experimental mechanical testing on 169 human samples to extract the elasticity properties of pathological breast tissues. The scale parameters were set to zero as these parameters are not incorporated in commercial finite element software packages like Comsol Multiphysics version 5.5. The thickness of the tissue samples was set to 1 cm while the tissue side length was 10 cm in this study. The Poisson’s ratio of all breast lesions was taken as 0.495 [
28] while the physical mass density of breast tissue was set to 945 kg/m
3 [
69,
70]. The patterns of in-plane wave propagation within the breast tissue, obtained using the finite element simulations, are also depicted in
Figure 4 for the first four in-plane modes. The visualisation was performed on the DCIS. Similar patterns were found for other breast lesion types, which are not included here for the sake of brevity.
Table 2 lists the mean square errors of training and testing datasets for the neural network of deep learning models studied in this paper. The outcome is indicated for three different deep learning architectures including (1) FL-3HDL(512,128,10), (2) FL-2HDL(128,10), and (3) FL-1HDL(10). Here, FL and HDL stand for the flatten layer and hidden dense layer, respectively. The integer preceding the ‘HDL’ indicates the quantity of hidden dense layers in the architecture of the neural network, and the corresponding numbers of neurons in each layer are denoted in parentheses in an ordered manner. For example, FL-2HDL(128,10) refers to a neural network with one flatten layer and two hidden dense layers, including 10 and 128 units within the second and first hidden layers, respectively. The mean square errors of these models are listed for various numbers of epochs from 1 to 500. The first model with one FL and three HDLs demonstrated the minimum mean square error, reaching around 4.71 × 10
−7 and 1.06 × 10
−6 after 500 epochs for the train and test datasets, respectively. In all cases, the mean square error of the deep learning model significantly decreased as the epoch number increased, indicating the importance of an appropriate epoch number in the accuracy of the machine learning modelling.
Table 3 compares the predicted dimensionless wave velocity within the early breast cancer of DCIS obtained using the deep learning model with those of the analytical approach (actual data) for different dimensionless (normalised) wave numbers. Scaling (normalisation) was conducted by dividing each feature by its maximum value to ensure that all features were in the same range.
The results of a number of different machine learning models are listed and compared in
Table 4. The machine learning models were used to predict the dimensionless wave velocity of in-plane elastic waves propagated within the breast lesions of women with high mammographic density. All models have been trained on 100 data points obtained using the higher-order nonlocal model. In addition to the neural network algorithm, Random Forest, Lasso, linear, Ridge, and ElasticNet regression models are presented in this table for comparison purposes. In this study, the neural network model consisted of one flatten layer and three hidden layers with 512, 128, and 10 neurons as well as one dense layer at the end with one neuron. The rectified linear unit (ReLU) function was used as the activation function to incorporate complexity in the form of nonlinearity and nonlocality into the model. The RandomForestRegressor from the open-source scikit-learn library was utilised with 100 estimators (decision trees). Each estimator made a prediction, and the final output was the average of these predictions for the regression problems. To adjust the machine learning models of regression and obtain the best hyperparameters, the GridSearchCV method of the sklearn.model_selection was used. The optimised values of the regularization parameter alpha for ElasticNet, Lasso, and Ridge were obtained as 0.005, 0.001, and 0.2, respectively. From this table, it is found that the neural network and Random Forest led to more accurate predictions compared to the other machine learning models. The results calculated using the neural network were slightly better than those of the Random Forest model (closer to the actual data) for various amounts of in-plane dimensionless wave numbers.
Figure 5 and
Figure 6 show the variation of phase velocity and frequency with respect to the in-plane wave number for the first and second modes, respectively. Healthy fibroglandular tissue, which is the main component of breast tissue, with high mammographic density and DCIS are compared in these figures. The first and second nonlocal parameters were set to 2 and 1 cm, respectively, and the strain gradient parameter was 0.5 cm. It was assumed that the thickness of the tissue was 1 cm. The phase velocity of the dense breast tissue was significantly lower than that of DCIS in both the first and second in-plane wave modes. In addition, the frequency of in-plane wave propagation within the DCIS was higher than that of the dense breast tissue due to the fact that the DCIS exhibited more rigidity under wave propagation and possessed a higher stiffness. In
Figure 5c,d, the influence of the length of the breast sample on the in-plane phase velocity curves is shown for human dense breast and DCIS, respectively. Various non-classical continuum models, together with the traditional elasticity model, were developed and used to plot these figures. In the higher-order nonlocal strain gradient model, the first and second nonlocal parameters were, respectively, set to 2 and 1 microns, and the strain gradient parameter was 0.5 microns. In the nonlocal model of the breast samples, the strain gradient effect was ignored while the first and second nonlocal parameters were 2 and 1 micron, respectively. By contrast, both the first and second nonlocal parameters were set to zero while the strain gradient parameter was taken as 0.5 micron in the strain gradient model. In the traditional classical model, the effects of all scale parameters were neglected, leading to a scale-free continuum model. The phase velocity estimated using the nonlocal model was lower than that of the classical model of breast tissues due to the stiffness-softening effect of stress nonlocality. By contrast, the strain gradient effect was associated with stiffness hardening and, hence, higher phase velocities. The nonlocal model incorporated stiffness softening at small scales while the strain gradient model could estimate structural stiffness hardening of the breast tissue. A refined combination of the models developed in this study simultaneously simulated both softening and hardening responses. At larger scales, as anticipated, all models exhibited convergence toward the classical model.
To study the influence of various scale parameters including first nonlocal, second nonlocal, and strain gradient on the in-plane wave propagation characteristics of breast lesions, phase velocity versus the in-plane wave number was plotted in
Figure 7a–c, respectively. The breast lesion was assumed to be of intermediate invasive ductal carcinoma (IDC) type, which accounts for around 75–80% of all invasive breast cancers [
71]. The elasticity and Poisson’s ratio of the intermediate IDC were taken as 19.99 kPa and 0.495, respectively [
28]. From
Figure 7a, it can be concluded that, when the first nonlocal parameter was increased from 0 to 2 cm, the phase velocity experienced a significant reduction. This change in phase velocity became higher by increasing the in-plane wave number. A slight reduction in the phase velocity of in-plane waves within the intermediate IDC by increasing the second nonlocal parameter was found (
Figure 7b). At small in-plane wave numbers (k < 30 1/m), the influence of the first-order stress nonlocality was negligible and the zeroth-order nonlocality was dominant.
Figure 7c depicts the variation of phase velocity with the wave number for different strain gradient parameters and for in-plane wave propagation. The phase velocity was very sensitive to the amount of strain gradient, especially at higher in-plane wave numbers. A small increase in the strain gradient parameter led to a sharp increase in phase velocity. Strain gradient effects were associated with stiffness hardening, leading to an increase in the frequency and phase velocity of in-plane waves within the human breast lesion.
To investigate and compare the wave propagation characteristics of various breast lesions, the frequency and phase velocity were plotted in
Figure 8 and
Figure 9 for nonlocal and strain gradient models, respectively. Poisson’s ratio of human breast tissue was considered to be 0.495. Elasticity moduli of adipose, fibroglandular, fibroadenoma, DCIS, low-grade, intermediate-grade, and high-grade IDC were 3.25 ± 0.91, 3.24 ± 0.61, 6.41 ± 2.86, 16.38 ± 1.55, 10.40 ± 2.60, 19.99 ± 4.2, and 42.52 ± 12.47, respectively [
28]. Furthermore, the elasticity moduli of the infiltrating lobular carcinoma (ILC), fibrocystic disease, invasive mucinous carcinoma (IMC), and fat necrosis were, respectively, 15.62 ± 2.64, 17.11 ± 7.35, 20.21 ± 4.2, and 4.45 ± 0.61. The in-plane wave number was set to 20 for this investigation. Based on the frequency and phase velocity, the differentiation between the common benign conditions such as the fibroadenoma and fat necrosis and cancerous breast lesions including IMC and IDC would be feasible in elastography biomedical imaging using in-plane elastic wave propagation. A remarkable difference was also observed between the in-plane wave response of healthy breast tissue and that of early breast cancer in the form of DCIS.
4. Discussion
From
Table 1, it can be concluded that there is an excellent match between the results of finite element simulation and those of the analytical model in all cases. The analytical model was developed based on the higher-order nonlocal theory, in which the gradients of strain were also incorporated. However, as the scale of the biological component increases, strain gradient and nonlocal effects become less important and can be ignored at a certain scale. In this case, the higher-order model reduces to the classical model, which can be validated using finite element simulations. The validation study was performed on three important breast tissues, including DCIS, healthy fibroglandular tissue, and fat necrosis. The DCIS and fibroglandular tissue are, respectively, taken into consideration as early breast cancers and a type of healthy tissue. In the mammography imaging of mammographically dense breasts, in which there is a significant amount of fibroglandular tissue compared to fat tissue, the DCIS is likely to be hidden by the fibroglandular tissue, making the diagnosis process challenging. In addition, fat necrosis was considered in this finite element study as a common type of benign condition.
Figure 4 shows that, at higher in-plane modes, the change in the tissue displacement occurs in smaller regions, making the deformation field more complex but more section-sensitive. Exciting higher vibration modes could be an effective way to improve the resolution and sensitivity of wave propagation-based detection tools like elastography imaging, as higher modes are associated with more ups and downs in the tissue, and, consequently, finer spatial details about the tissue’s biomechanical features.
Table 2 indicates the importance of an appropriate architecture design for the deep learning model. When the number of hidden layers is increased from one to three, the MSE of both test and train data significantly decreases. This significant reduction in the MSE of the model is not only related to the number of hidden layers but also is associated with the number of neurons in the architecture. The third model of deep learning, which is called FL-1HDL(10) in
Table 2, contains 10 neurons, while the first model with the name FL-3HDL(512,128,10) incorporates 650 neurons.
Table 3 further confirms this finding by comparing the results of the deep learning model with those of actual data obtained by using the higher-order nonlocal model of elasticity. It is observed that a deep learning algorithm with one flatten layer and two dense layers containing around several hundred neurons can highly accurately predict the results of the nonlocal elasticity model without having any elasticity theory knowledge.
The rigid-tumour–soft-cell paradox is associated with stiffness hardening at the tissue scale level while the stiffness of a clump of cells experiences a softening behaviour. The classical model of elasticity can only predict stiffness hardening and lacks the simulation of stiffness softening at cell scales. The nonlocal elasticity model with higher-order scale and strain parameters can predict both stiffness hardening and softening depending on the scale of the biological component or tissue. At large scales, the strain gradient effect is dominant, which is highly linked with stiffness hardening, as can be concluded from
Figure 5,
Figure 6 and
Figure 7. However, at certain cell levels, the stress nonlocalities of zeroth and first orders account for the stiffness softening, as observed in
Figure 7a,b. Specifically, in terms of wave propagation characteristics, stress nonlocalities lead to a reduction in stiffness, and this results in lower frequencies and phase velocities, which is supported by the data in
Figure 8 and
Figure 9. By contrast, strain gradients are linked to a significant enhancement in the total stiffness, leading to higher in-plane wave frequencies and phase velocities. Another important finding is that, when the wave mode number increases, the difference between the malignant and healthy breast tissue significantly increases at both tissue and cellular levels, making the frequency-based diagnosis process easier and more accurate. Moreover, stress nonlocalities become more prominent at higher in-plane modes since the interaction between different sections of the breast tissue is greater when the in-plane wave number increases. Based on the comparison study made in this research, a significant difference between the in-plane wave propagation characteristics of healthy fibroglandular/adipose tissues and that of early breast cancer (DCIS) was observed, which makes wave-propagation-based imaging tools like elastography promising and feasible for the early detection of breast cancer. Unlike X-ray mammography, there is no masking effect when the density (i.e., the relative amount of fibroglandular tissue) of the human breast increases, making in-plane wave propagation-based elastography imaging a promising complementary imaging technology for women with dense breasts.
5. Conclusions
In this study, a higher-order model of nonlocal elasticity in conjunction with a deep learning algorithm was developed and used to accurately obtain the in-plane wave propagation characteristics of human breast tissues. Eleven different healthy, benign, and pathological conditions including adipose, fibroglandular, DCIS, ILC, IDC (intermediate, low- and high-grade), IMC, fat necrosis, and fibrocystic and fibroadenoma were taken into account. Using Hamilton’s law, the differential equations of the in-plane motion of breast tissue were derived. Analytical solutions were presented for the frequency and phase velocity of the elastic waves propagated within the thin breast tissue sections. To show the validity of the model, three finite element studies were conducted on DCIS as the very early form of breast cancer, using fibroglandular tissue as the healthy mammographically dense breast and fat necrosis as a common benign condition. The results of the finite element studies were in excellent agreement with those of higher-order nonlocal elasticity modelling. Furthermore, an artificial deep neural network was designed, created, and adjusted for the in-plane wave propagation analysis of breast lesions. The deep learning model can discover the mathematical formulas that govern the in-plane wave propagation response of breast lesions and precisely predict the phase velocity with training on only a few hundred samples. For the deep learning model with FL-3HDL(512,128,10) architecture, the mean square error of training and testing datasets were, respectively, reduced to 4.7 × 10−7 and 1.06 × 10−6 after 500 epochs. The results of the deep learning model were also compared to those of other machine learning models, including linear, Ridge, Lasso, Random Forest, and ElasticNet regression models. An excellent match was found between the predictions made by the neural network algorithm and that of the Random Forest model. The phase velocity of in-plane waves within the DCIS was, on average, more than 2.5 times higher than that of healthy fibroglandular tissue. Moreover, a significant difference (p-value < 0.01) was found between breast cancers such as high-grade IDC and IMC and benign conditions like fat necrosis and fibroadenoma. Both zeroth- and first-order stress nonlocalities were associated with a reduction in the phase velocity of in-plane waves since they were linked to tissue stiffness softening. By contrast, the strain gradient parameter enhanced the in-plane phase velocity of breast tissue as it was correlated with stiffness hardening. The combination of nonlocal and strain gradient parameters in the higher-order elasticity model allowed for the concurrent incorporation of stiffness hardening and softening, solving the rigid-tumour–soft-cell paradox of breast biomechanics.