Flexural Eigenfrequency Analysis of Healthy and Pathological Tissues Using Machine Learning and Nonlocal Viscoelasticity

Abstract

Biomechanical characteristics can be used to assist the early detection of many diseases, including breast cancer, thyroid nodules, prostate cancer, liver fibrosis, ovarian diseases, and tendon disorders. In this paper, a scale-dependent viscoelastic model is developed to assess the biomechanical behaviour of biological tissues subject to flexural waves. The nonlocal strain gradient theory, in conjunction with machine learning techniques such as extreme gradient boosting, k-nearest neighbours, support vector machines, and random forest, is utilised to develop a computational platform for biomechanical analysis. The coupled governing differential equations are derived using Hamilton’s law. Transverse wave analysis is conducted to investigate different normal and pathological human conditions including ovarian cancer, breast cancer, and ovarian fibrosis. Viscoelastic, strain gradient, and nonlocal effects are used to describe the impact of fluid content, stiffness hardening caused by the gradients of strain components, and stiffness softening associated with the nonlocality of stress components within the biological tissues and cells. The integration of the scale-dependent biomechanical continuum model with machine learning facilitates the adoption of the developed model in practical applications by allowing for learning from clinical data, alongside the intrinsic mechanical laws that govern biomechanical responses.

1. Introduction

The mechanical properties of cells, as well as the mechanical characteristics of their microenvironment, affect cell behaviour, differentiation, and cell fate [1]. Understanding the interaction between mechanical signals and biological properties plays a vital role in the appropriate design of tissue engineering scaffolds and the advancement of regenerative medicine [2,3]. It has been revealed that the regenerative characteristics of mesenchymal stem cells can be enhanced by means of mechanical conditioning [4]. Moreover, due to the importance of mechanical stimuli in the tumour microenvironment, understanding the mechanobiological characteristics of cancer cells would enhance more targeted cancer treatments and enable strategies for cancer management [5]. In the clinic, the mechanical response of human tissue under wave propagation induced by external transducers has recently been utilised to assist in the detection of liver fibrosis [6], cardiovascular risk [7], breast cancer [8], and ovarian diseases [9].

To assess the biomechanical and mechanobiological properties of both biological tissues and individual cells, mathematical models are of high importance, as they bridge the gap between experimental observations and theoretical predictions [10]. These models provide valuable insights into the behaviour and responses of complex biological systems subject to mechanical stimuli. Particularly important to disease detection, continuum theories such as poroelasticity, elasticity, and viscoelasticity enable us to develop accurate computational and analytical platforms to extract intrinsic biomechanical features, which are essential for early disease diagnosis where tissue stiffness hardening or softening is involved [11,12]. In addition to continuum mechanics-based models, molecular dynamics approaches have been widely used to assess the dynamic behaviour of sub-cellular components such as proteins, lipid membranes, organelles, and the cytoskeleton [13]. The application of an appropriate theoretical framework is highly dependent on the scale at which the mechanical properties are examined. At small-scale levels, such as nanoscale and microscale, molecular dynamics is used [14], while at the tissue level, finite element methods [15] and mathematical models of continuum mechanics [16] are applicable. In molecular dynamics, the motions of individual molecules over time are simulated and monitored by simulating their interactions with other molecules [13]. This simulation allows us to understand the influence of various factors and signals on the deformation and mechanical behaviour of proteins, biomolecules, and organelles. In contrast, at larger scales, including tissue and organ levels, classical continuum theories are used to develop analytical mathematical models to investigate biomechanical properties such as stiffness, vibrational frequency, and wave propagation response. Finite element techniques, in which the organ is divided into small, interconnected pieces, are commonly used where there is complexity in the analysis due to geometrical and material nonlinearities, as well as structural irregularities [17].

Size-dependent models of nonlocal continuum mechanics have shown great promise for understanding mechanical behaviour at intermediate small-scale levels between the molecular level and large-scale levels, and in bridging the gap between molecular dynamics and traditional continuum theories [18]. In advanced nonlocal theories of elasticity, extra mathematical terms are introduced, allowing for the incorporation of a broader spectrum of phenomena at small-scale levels, which are omitted in classical elasticity theories [19]. These phenomena include, but are not limited to, deformations due to higher gradients of strain components [20] and stiffness hardening associated with couple stresses [21], as well as stiffness softening induced by nonlocal mechanical stresses [22].

Gao and Lei [23] developed a size-dependent Timoshenko beam model based on the nonlocal elasticity theory to investigate the buckling behaviour of a single protein microtubule in a viscoelastic surrounding cytoplasm. Furthermore, Heireche et al. [24] introduced a modified beam model to analyse the vibrational characteristics of protein microtubules in living cells by incorporating the influences of nonlocal stress, rotary inertia, and shear deformation. Civalek and Demir [25] developed a scale-dependent finite element technique to assess the mechanical response of protein microtubules embedded in an elastic matrix. Moreover, Akgöz and Civalek [26] presented a higher-order shear deformation model of microtubules based on the modified strain gradient theory. They found that the influence of shear deformation was more pronounced when the aspect ratio of the microtubule was smaller, and nonlocal effects were higher for smaller lengths. Furthermore, the torsional and axial frequency responses of microtubules have been investigated according to a nonlocal continuum rod model [27]. In addition to classical nonlocal continuum mechanics [28], other scale-dependent models, including those based on nonlocal strain gradient theory [29,30], nonlocal couple stress theory [31], and surface elasticity theory [32], as well as nonlocal integral models [33], have been utilised to explore the biomechanical behaviour of biological components at small-scale levels. Moreover, a nonlocal thermoelasticity approach [34], a higher-order nonlocal theory [35], a nonlinear nonlocal elasticity model [36], and a coarse-grained differential-tension technique [37] have recently been developed for analysing the biomechanics of different biological tissues, indicating the validity and great promise of nonlocal scale-dependent models in biomechanical modelling.

In the present paper, we aim to develop a scale-dependent viscoelastic model for flexural wave propagation within biological cells and tissues. The proposed modified continuum model is capable of describing flexural wave characteristics at different scales, as both stiffness hardening and softening are considered within the framework of the nonlocal strain gradient model. The strain gradient parameter is responsible for the stiffness hardening response under external mechanical stimuli at larger scales, while the nonlocal stress parameter accounts for stiffness softening at smaller scales. The Kelvin–Voigt model of structural viscoelasticity is employed to incorporate the influence of the fluid content of biological tissue. Differential equations of wave dispersion are derived based on Hamilton’s law. To overcome the limitations of the nonlocal continuum mechanics approach caused by fundamental assumptions, a machine learning algorithm is developed and trained on the flexural wave characteristics of ovarian tissue with different fibrosis conditions as well as ovarian cancer and breast cancer. This work introduces the first viscoelastic higher-order nonlocal model to investigate flexural wave propagation in ovarian tissue. The machine learning algorithm demonstrates great potential in accurately identifying the eigenfrequencies and damping ratios of ovarian and breast tissues under flexural waves, even with a limited number of samples. Using the viscoelastic model of nonlocal elasticity, significant differences are identified between fibrous ovarian tissue and normal tissue, as well as between the ovarian cancerous and healthy condition, demonstrating that the eigenfrequency parameter and damping ratio can serve as indicators of ovarian diseases.

2. Scale-Dependent Continuum Approach of Nonlocal Viscoelasticity

Assuming in-plane displacements are small, the strain components of a biological sample (${\epsilon }_{ij}$) can be expressed in terms of the transverse displacement (w) as

$$\begin{gathered} {\epsilon }_{xx}(x,y,z,t)=-z{\displaystyle \frac{{\partial }^{2}w(x,y,t)}{\partial x^2}}, \\ {\epsilon }_{yy}(x,y,z,t)=-z{\displaystyle \frac{{\partial }^{2}w(x,y,t)}{\partial y^2}}, \\ {\epsilon }_{xy}(x,y,z,t)=-z{\displaystyle \frac{{\partial }^{2}w(x,y,t)}{\partial x\partial y}}, \end{gathered}$$

(1)

where x, y, and z are the components of the Cartesian coordinate system and t is time. Based on nonlocal strain gradient continuum mechanics [20], the components of the couple stress resultants ($M_{ij}$) are obtained:

$$\begin{array}{l}\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)M_{xx}=\\ -D_{11}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\right)\right]\\ -D_{12}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)\right],\\ \left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)M_{yy}=\\ -D_{12}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\right)\right]\\ -D_{22}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)\right],\\ \left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)M_{xy}=-2D_{33}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial x\partial y}}-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial x\partial y}}\right)\right],\end{array}$$

(2)

where

$$\begin{gathered} {\nabla }^{2}f={\displaystyle \frac{d^{2}f}{d{x}^2}}+{\displaystyle \frac{d^{2}f}{d{y}^2}}+{\displaystyle \frac{d^{2}f}{d{z}^2}}, M_{ij}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}z{\sigma }_{ij}dz}, \\ \begin{array}{l}\hspace{1em} D_{11}={\displaystyle \frac{E_{11}{h}^3}{12(1-{\nu }_{12}{\nu }_{21})}}, D_{22}={\displaystyle \frac{E_{22}{h}^3}{12(1-{\nu }_{12}{\nu }_{21})}}, \\ \hspace{1em} D_{12}=D_{21}={\displaystyle \frac{{\nu }_{12}{E}_{22}{h}^3}{12(1-{\nu }_{12}{\nu }_{21})}}, D_{33}={\displaystyle \frac{G_{12}{h}^3}{12}},\end{array} \end{gathered}$$

(3)

in which ${\nabla }^{2}f$ is the Laplacian operator of an arbitrary given function f. The thickness of the rectangular-shaped biological sample is indicated by h. The elasticity modulus, shear modulus, and Poisson’s ratio of the sample are denoted by E_ij, G_ij and v_ij, respectively. ${\lambda }_0^{nls}$ and ${\lambda }_1^{nls}$ represent the scale parameters of the first and second stress nonlocalities, which are defined by ${\lambda }_0^{nls}={\left({\kappa }_0{\mathcal{l}}_{cell}\right)}^2$ and ${\lambda }_1^{nls}={\left({\kappa }_1{\mathcal{l}}_{cell}\right)}^2$, respectively [20,30]. In addition, ${\mathsf{\Theta }}_{stg}^2$ is the strain gradient parameter, which is also indicated by ${\mathcal{l}}_s^2$ (i.e., ${\mathsf{\Theta }}_{stg}^2={\mathcal{l}}_s^2$). ${\kappa }_0$ and ${\kappa }_1$ are two calibration parameters used to better mimic experimental data. ${\mathcal{l}}_{cell}$ is a cellular-level geometrical parameter used to incorporate the effect of the inner structural organisation of the biological sample. Based on the Kelvin–Voigt model of viscoelasticity, Young’s modulus E_ij is replaced by $E_{ij}\left[(*)+g\left(\partial /\partial t\right)(*)\right]$, and consequently, the stiffness component D_ij is replaced by $D_{ij}\left[(*)+g\left(\partial /\partial t\right)(*)\right]$, where g is the viscoelastic damping coefficient [38]. The final couple stress resultants of the scale-dependent continuum approach with viscoelasticity are given by

$$\begin{array}{l}\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)M_{xx}=\\ -D_{11}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial x^2}}\right)-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial x^2}}\right)\right]\\ -D_{12}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial y^2}}\right)-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}+g{\displaystyle \frac{{\partial }^{2}w}{\partial t\partial y^2}}\right)\right],\\ \left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)M_{yy}=\\ -D_{12}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial x^2}}\right)-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial x^2}}\right)\right]\\ -D_{22}\left[\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial y^2}}\right)-{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial y^2}}\right)\right],\\ \left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)M_{xy}=-2D_{33}\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left({\displaystyle \frac{{\partial }^{2}w}{\partial x\partial y}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial x\partial y}}\right)\\ +2D_{33}{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left({\displaystyle \frac{{\partial }^{2}w}{\partial x\partial y}}+g{\displaystyle \frac{{\partial }^{3}w}{\partial t\partial x\partial y}}\right).\end{array}$$

(4)

Using Hamilton’s law, the equation of the time-dependent transverse deformation of the biological sample can be written as

$$\begin{array}{l}{\displaystyle \frac{{\partial }^2{M}_{xx}}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{M}_{yy}}{\partial y^2}}+2{\displaystyle \frac{{\partial }^2{M}_{xy}}{\partial x\partial y}}+f_z\\ =m_0{\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}-m_2\left({\displaystyle \frac{{\partial }^{4}w}{\partial x^2\partial t^2}}+{\displaystyle \frac{{\partial }^{4}w}{\partial y^2\partial t^2}}\right),\end{array}$$

(5)

where m₀ and m₂ are the mass inertia terms and are related to the mass density of the biological sample (${\rho }_{bio}$) by the following relationships:

$$m_0={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{\rho }_{bio}dz}, m_2={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}{z}^2{\rho }_{bio}dz}.$$

(6)

It is assumed that the sample is embedded in a viscoelastic medium with normal, shear, and structural viscoelasticity. The external transverse load per unit area (f_z) due to the surrounding viscoelastic medium is

$$f_z=-k_{w}w+k_s\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)-c_f{\displaystyle \frac{\partial w}{\partial t}},$$

(7)

Here, k_w, k_s, and c_f are the Winkler constant, shear constant, and viscoelastic coefficient of the surrounding medium, respectively. Substituting the couple stress components from Equation (4) into Equation (5), and using Equation (7), one obtains

$$\begin{array}{l}-\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left[D_{11}\left({\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}+g{\displaystyle \frac{{\partial }^{5}w}{\partial t\partial x^4}}\right)+D_{22}\left({\displaystyle \frac{{\partial }^{4}w}{\partial y^4}}+g{\displaystyle \frac{{\partial }^{5}w}{\partial t\partial y^4}}\right)\right.\\ \left.+2\left(D_{12}+2D_{33}\right)\left({\displaystyle \frac{{\partial }^{4}w}{\partial x^2\partial y^2}}+g{\displaystyle \frac{{\partial }^{5}w}{\partial t\partial x^2\partial y^2}}\right)\right]\\ +{\mathsf{\Theta }}_{stg}^2\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\nabla }^2\left[D_{11}\left({\displaystyle \frac{{\partial }^{4}w}{\partial x^4}}+g{\displaystyle \frac{{\partial }^{5}w}{\partial t\partial x^4}}\right)+D_{22}\left({\displaystyle \frac{{\partial }^{4}w}{\partial y^4}}+g{\displaystyle \frac{{\partial }^{5}w}{\partial t\partial y^4}}\right)\right.\\ \left.+2\left(D_{12}+2D_{33}\right)\left({\displaystyle \frac{{\partial }^{4}w}{\partial x^2\partial y^2}}+g{\displaystyle \frac{{\partial }^{5}w}{\partial t\partial x^2\partial y^2}}\right)\right]\\ +\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)\left[-k_{w}w+k_s\left({\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}w}{\partial y^2}}\right)-c_f{\displaystyle \frac{\partial w}{\partial t}}\right]\\ =m_0\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right){\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}\\ -m_2\left(1-{\lambda }_1^{nls}{\nabla }^2\right)\left(1-{\lambda }_0^{nls}{\nabla }^2\right)\left({\displaystyle \frac{{\partial }^{4}w}{\partial x^2\partial t^2}}+{\displaystyle \frac{{\partial }^{4}w}{\partial y^2\partial t^2}}\right).\end{array}$$

(8)

The transverse deflection of the biological sample induced by the propagation of flexural waves can be expressed as

$$w\left(x,y,t\right)=W\exp \left[i\left(k_{x}x+k_{y}y-\omega t\right)\right],$$

(9)

where W is a constant associated with the transverse deflection amplitude. k_x and k_y are wave numbers in the x and y directions, respectively. $\omega$ represents the dimensional frequency parameter. Substituting Equation (9) into the time-dependent equation of deformation leads to

$$\begin{array}{l}{\gamma }_{12}^{(-)}\left[D_{11}{\left(k_x\right)}^4+D_{22}{\left(k_y\right)}^4+2\left(D_{12}+2D_{33}\right){\left(k_x\right)}^2{\left(k_y\right)}^2\right]\\ +{\gamma }_{10}^{(\times )}\left\{k_w+k_s\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]\right\}\\ -i\omega \left\{c_f{\gamma }_{10}^{(\times )}+g{\gamma }_{12}^{(-)}\left[D_{11}{\left(k_x\right)}^4+D_{22}{\left(k_y\right)}^4+2\left(D_{12}+2D_{33}\right){\left(k_x\right)}^2{\left(k_y\right)}^2\right]\right\}\\ ={\gamma }_{10}^{(\times )}\left\{m_0+m_2\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]\right\}{\omega }^2.\end{array}$$

(10)

where

$$\begin{array}{l}{\gamma }_{12}^{(-)}=1+{\lambda }_1^{nls}\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]+{\mathsf{\Theta }}_{stg}^2\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]\left(1+{\lambda }_0^{nls}\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]\right),\\ {\gamma }_{10}^{(\times )}=\left(1+{\lambda }_1^{nls}\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]\right)\left(1+{\lambda }_0^{nls}\left[{\left(k_x\right)}^2+{\left(k_y\right)}^2\right]\right).\end{array}$$

(11)

For generalisation and simplification, a set of dimensionless geometrical and physical parameters is introduced as

$$\begin{array}{l}{\overline{k}}_w={\displaystyle \frac{k_w{L}^4}{D_{11}}}, {\overline{k}}_s={\displaystyle \frac{k_s{L}^2}{D_{11}}}, {k^{\prime }}_x=k_x{L}_x, {k^{\prime }}_y=k_y{L}_y,\\ {\chi }_0={\displaystyle \frac{e_{0}a}{L_x}}, {\chi }_1={\displaystyle \frac{e_{1}a}{L_x}}, {\chi }_s={\displaystyle \frac{{\mathcal{l}}_s}{L_x}}, r={\displaystyle \frac{L_x}{L_y}}, L_x=L, \\ \overline{\omega }=\omega L^2\sqrt{{\displaystyle \frac{{\rho }_{bio}h}{D_{11}}}}, Z_{ij}={\displaystyle \frac{D_{ij}}{D_{11}}}, {\overline{c}}_f={\displaystyle \frac{c_f{L}^2}{\sqrt{{\rho }_{bio}h{D}_{11}}}}, \overline{g}={\displaystyle \frac{g}{L^2}}\sqrt{{\displaystyle \frac{D_{11}}{{\rho }_{bio}h}}},\\ \overline{D}={\left({k^{\prime }}_x\right)}^4+2r^2\left(Z_{12}+2Z_{33}\right){\left({k^{\prime }}_x\right)}^2{\left({k^{\prime }}_y\right)}^2+r^4{Z}_{22}{\left({k^{\prime }}_y\right)}^4,\\ {\overline{k}}_{ws}={\overline{k}}_w+{\overline{k}}_s\left[{\left({k^{\prime }}_x\right)}^2+r^2{\left({k^{\prime }}_y\right)}^2\right], \overline{m}=1+{\displaystyle \frac{h^2}{12L^2}}\left[{\left({k^{\prime }}_x\right)}^2+r^2{\left({k^{\prime }}_y\right)}^2\right],\end{array}$$

(12)

where $\overline{\omega }$ is the dimensionless frequency parameter. In this analysis, viscoelastic effects associated with the fluid content of both the biological sample and the surrounding medium are incorporated. In viscoelasticity, the frequency parameter is divided into real and imaginary parts as

$$i\overline{\omega }=\mathsf{\Omega }\left(\zeta \pm i\sqrt{1-{\zeta }^2}\right),$$

(13)

where $\zeta$ is the general damping ratio and $\mathsf{\Omega }$ is the undamped frequency parameter of the biological system.

Table 1. Limitations of the current modelling approach for the estimation of biomechanical characteristics compared to available models.

Model	Model’s Biomechanical Parameters	Limitations	Scale Effects
Classical elasticity	Elastic constants, Poisson’s ratio	Linearity, homogeneity, homogeneity, boundary conditions	Not incorporated
Classical poroelasticity	Elastic moduli, Poisson’s ratio, fluid-related parameters	Simplified fluid interaction, homogeneity, coupling complexity	Not incorporated
Local viscoelasticity	Elastic moduli, Poisson’s ratio, viscoelastic damping coefficient	Linear assumption, homogeneity, boundary conditions	Not incorporated
Nonlocal elasticity	Elastic moduli, Poisson’s ratio, stress nonlocality	Parameter identification, homogeneity, viscoelastic effects	Zeroth-order nonlocal effect
Present model	Elastic moduli, Poisson’s ratio, viscoelastic damping, strain gradient, stress nonlocality	Computational complexity, parameter identification, homogeneity	Three different scale effects

lists the limitations of the proposed viscoelastic model of nonlocal strain gradient theory compared to other continuum models. The classical elasticity model is commonly used in its linear form, which is restricted to small deformations. Furthermore, classical elasticity models are not able to describe fluid-related effects such as viscoelastic damping and fluid–solid interactions. The boundary conditions of biological tissues in practical situations are complex, making their implementation challenging and computationally expensive. Classical poroelasticity theory assumes a simplified fluid flow model within a porous structure, which might not be able to completely model complex fluid–matrix interactions in biological tissues. The classical poroelasticity theory lacks scale effects, restricting its capability for detecting tumours to small-scale levels. The application of local viscoelasticity models in the detection of abnormalities is limited due to homogeneity and scale invariance assumptions. The nonlocal elasticity theory incorporates a scale parameter that is associated with stiffness softening caused by stress nonlocality. The present model of nonlocal strain gradient viscoelasticity contains three scale parameters, allowing for the comprehensive incorporation of scale effects. However, the computational cost of the present model is higher than those available in the literature. Furthermore, the precise determination of these scale parameters from experimental data can be challenging and costly in healthcare applications.

3. Finite Element Approach

COMSOL Multiphysics Simulation software version 5.5 was employed to build a finite element model (FEM) to simulate the mechanical behaviour of biological tissues under flexural wave propagation at large-scale levels and to generate more data for training a machine learning model. The three-dimensional shell interface of the structural mechanics module was utilised to create and analyse the tissue. Shells are thin structures with flexural stiffness that can be either flat or curved. In this analysis, the shell structure was assumed to be flat and of a rectangular shape. Geometrical nonlinearity was assumed to be negligible due to the application of light flexural waves in healthcare. To solve the FEM-based governing of differential equations, a multifrontal massively parallel sparse direct (MUMPS) technique was implemented. The sequence type to create finite elements was set to the physics-controlled mesh with an extra fine element size. The four edges of the section were assumed to be simply supported. The initial displacements and velocity components were set to zero. The length, width, and thickness were taken as 30 cm, 30 cm, and 1 cm, respectively. The effects of the shear modulus on the flexural wave response were negligible, as the length-to-thickness ratio was considerably high.

4. Machine Learning Approach

In this section, several machine learning algorithms are developed for analysing the biomechanical flexural wave characteristics of biological tissues. The advantages of machine learning approaches, with respect to other techniques such as scale-dependent continuum mechanics, molecular dynamics, and finite element methods, include, but are not limited to, scalability [39], computational efficiency [40], (bio-) complexity handling [41], cyber threat management [42], and the ability to analyse large datasets in real time [43]. Each mathematical model developed based on continuum mechanics comes with its own assumptions that restrict its application in practical situations. Moreover, simulation techniques such as the finite element method (suitable for large-scale modelling) and molecular dynamics (suitable for small-scale modelling) are particularly expensive, restricting their applications in the clinic. By contrast, the use of a properly trained machine learning model in healthcare for the extraction of biomechanical features of biological samples at different scales is rapid and scalable to large datasets [44,45]. In addition, machine learning excels in the incorporation of complexities such as nonlinear responses, non-uniform loading conditions, and imperfections. These capabilities and benefits make machine learning a desirable approach to study the flexural wave characteristics of biological samples.

Figure 1 provides a brief overview of the proposed machine learning (ML) methodology for studying the biomechanical characteristics of biological samples. The ML tool utilises a hybrid dataset for training, including scale-dependent data derived from the nonlocal strain gradient model of viscoelasticity, experimental data, and flexural wave characteristics based on the FEM. The biomechanical properties of various biological samples, including both control samples and pathological conditions, are used during the training process to fit a regression-based ML algorithm. A prediction is made to simulate the frequencies of various samples under flexural waves. To verify the methodology’s accuracy and efficacy, mean squared errors (MSE) are obtained by comparing the predicted frequencies with those of the test dataset. Train/test splitting is performed considering 70% of the data for training and 30% for testing.

Different ML regressors are developed and tested, including Ridge, ElasticNet, Lasso, random forest [46], and extreme gradient boosting (XGBoost), as well as the regressors of the support vector machine (SVM) technique [47] and k-nearest neighbours (KNN). The present integrated methodology harnesses the advantages of ML regressors in conjunction with fundamental scale-dependent principles of viscoelasticity and experimental observation, providing a comprehensive framework for understanding tissue behaviour under transverse waves. The KNN, SVM, and random forest regressors, together with Ridge, Lasso, and ElasticNet, are imported from scikit-learn library 1.2.2 [48]. To normalise the features, they are divided by their corresponding maximum value. The target is set as the frequency parameter of the biological tissue under flexural wave propagation. A biomechanical dataset with 201 biological samples for each case study is employed. In healthcare and biomedical applications, providing a large dataset with different types of malignancy and healthy conditions for training is often a challenge. Therefore, data augmentation techniques and machine learning models with the capability of handling small datasets effectively are of high importance. The integration of the ML model with the FEM and nonlocal strain gradient model of viscoelasticity would further enhance the diversity of the dataset and help prevent overfitting.

5. Results

5.1. Validation Study

FEM simulations were conducted to verify the results of the scale-dependent viscoelastic model considering three different case studies: (1) human ovarian cancer at the cellular level, (2) human breast cancer at the cellular level, and (3) mouse ovarian fibrosis associated with aging. The mechanical properties of OVCAR-3 (ovarian cancer cells), HO-8910 (ovarian cancer cells), and HOSEpiC (ovarian healthy cells) were extracted from the first FEM analysis in Ref. [49]. Two types of human breast cells, including MCF-7 as cancer cells and MCF-10A as healthy cells with a known Young’s modulus [50], were considered in the FEM simulations of the second case study. Furthermore, the eigenfrequencies of the flexural waves within different ovarian mouse tissues were investigated with a particular focus on ovarian fibrosis with aging. The following ovarian mouse tissues were simulated: old, young, and old with collagenase treatment. The Young’s moduli of these mouse tissues were experimentally measured using instrumental indentation [51]. Nonlocal and strain gradient effects were not incorporated, since COMSOL Multiphysics Simulation software version 5.5 did not encompass the nonlocal strain gradient theory.

Table 2. Comparison study: FEM versus the present scale-dependent viscoelastic model. Ovarian cancer cells of type OVCAR-3 and HO-8910, breast cancer cells of type MCF-7, as well as healthy ovarian (HOSEpiC) and breast (MCF-10A) cells are taken into account. Furthermore, a comparison study is conducted between the eigenfrequencies of different ovarian mouse tissues.

Type	Young’s Modulus (Pa)	Dimensionless Flexural Wave Number	FEM	Present Model	Percentage Error (%)
OVCAR-3 (ovarian cancer cells)	1195.72	1	0.1297	0.1305	0.6168
	1195.72	2	0.5104	0.5218	2.2335
HO-8910 (ovarian cancer cells)	996.27	1	0.1184	0.1191	0.5912
	996.27	2	0.4659	0.4763	2.2322
HOSEpiC (ovarian healthy cells)	2160.94	1	0.1744	0.1754	0.5734
	2160.94	2	0.6861	0.7015	2.2446
MCF-7 (breast cancer cells)	487.44	1	0.0828	0.0833	0.6039
	487.44	2	0.3259	0.3332	2.2400
MCF-10A (breast healthy cells)	1231.07	1	0.1316	0.1324	2.2400
	1231.07	2	0.5179	0.5295	2.2398
Young mouse ovary	1980	1	0.1669	0.1679	0.5992
		2	0.6567	0.6715	2.2537
Old mouse ovary	4360	1	0.2477	0.2491	0.5652
		2	0.9746	0.9964	2.2368
Old mouse ovary with collagenase treatment	2280	1	0.1791	0.1801	0.5583
		2	0.7047	0.7205	2.2421

lists the eigenfrequencies obtained by the present scale-dependent viscoelastic model and those estimated by the FEM. The percentage error was defined as 100 × (FEM − NSGT)/FEM, where NSGT is an abbreviation for nonlocal strain gradient theory. A close match was found between the FEM and the NSGT for the eigenfrequencies of vibrations induced by the propagation of flexural waves in all three case studies.

The results of the present viscoelastic model are compared with those reported in the literature on the scale-dependent eigenfrequency response of small-scale structures [38] in Figure 2. The present scale-dependent viscoelastic model contains three different scale parameters, including two nonlocal parameters and one strain gradient parameter. When the strain gradient effects are ignored and the scale influence of the second nonlocal parameter is neglected, the present model is reduced to the classical nonlocal viscoelastic model [38]. To perform a reasonable comparison, the dimensionless material and geometrical properties available in the model are set the same as in Ref. [38]. The dimensionless Winkler coefficient of the surrounding medium is taken as 100, while its shear coefficient is set to zero. The validation study is conducted for two different values of the nondimensional damping coefficients (NDCs) of the surrounding medium (i.e., NDC = 20 and NDC = 28). It is found that the results of the present viscoelastic model are in excellent agreement with those reported in the open literature on small-scale structures.

5.2. Ovarian Cancer

The first case study involves the biomechanical analysis of ovarian cancer. Recently, it has been shown that the biomechanical features of ovarian cells can be used as indicators of ovarian cancer state and progression [49]. Three different biological conditions are taken into account: (1) OVCAR-3 human ovarian cancer cells, (2) HO-8910 ovarian cancer cells, and (3) human ovarian surface epithelial cells (HOSEpiC) as the healthy (control) condition. The viscoelastic and Young’s moduli of these cells are extracted from the available data reported in the literature [49], which had been experimentally obtained using scanning force microscopy. The average elasticity moduli of the OVCAR-3, HO-8910, and HOSEpiC cells are 1195.72, 996.27, and 2160.94 Pa, respectively. The structural viscoelasticity coefficient of the ovarian cancer cells is 11.38 Pa·s for OVCAR-3 and 10.70 Pa·s for HO-8910, while the healthy ovarian cells (HOSEpiC) have a viscoelasticity coefficient of 30.00 Pa·s. The Poisson’s ratio and mass density are, respectively, set to 0.495 and 945 kg/m³ [15,52,53]. It is assumed that the biomechanical properties of the ovarian cells are isotropic. The ratio of the side length of the sample to the sample thickness is set to 10 while the ratio of the sample length to its width is equal to one. The wave numbers in the x and y directions are assumed to be the same. For visualisation, the non-dimensional wave numbers in the x and y directions are obtained by $k_x{L}_x/\pi$ and $k_y{L}_y/\pi$, respectively. It is assumed that the ovarian cell samples are embedded in a viscoelastic medium made of engineered dextran-based hydrogels with an elasticity modulus of 16 kPa and a viscoelastic coefficient of 4 kPa·s [54]. The ratio of the thickness of the viscoelastic surrounding medium to the thickness of the sample is taken as 0.1. An intermediate value of 0.05 is used for all dimensionless scale parameters in the nonlocal strain gradient model.

Figure 3 indicates the imaginary part of the eigenfrequency for the HOSEpiC, HO-8910, and OVCAR-3 cells versus the flexural wave number. The imaginary part of the eigenfrequency gradually reduces as the flexural wave number increases until it reaches its critical point, in which the imaginary part is zero. The healthy ovarian cells are the first sample that reach their critical point, at around the dimensional wave number 7.4. By contrast, the critical wave number where the imaginary part of the eigenfrequency of ovarian cancer cells becomes zero is around 10. Figure 4 depicts the real part of the eigenfrequency of different ovarian samples versus the dimensionless wave number. It is observed that, at larger wave numbers, the difference between the frequency response of ovarian cancer and healthy cells is more pronounced. Before the critical point, there is only one trajectory for the real part of the eigenfrequency of flexural waves within ovarian cells, while it divides into two branches after the critical point. In all pathological cases, the lower branch is marked with an asterisk. The structural damping ratio and undamped frequency parameter are also plotted in Figure 5 and Figure 6 against the flexural wave number for the ovarian healthy cells (HOSEpiC) and ovarian cancer cells (HO-8910 and OVCAR-3), respectively. The damping ratios of both healthy and cancer cells increase with the increase in wave number. Furthermore, the difference between the damping ratios of healthy cells and ovarian cancer cells becomes greater by increasing the wave number from 0 to 10. At the dimensionless wave number of 10, the damping ratio of HOSEpiC is about 2.5 times higher than that of the ovarian cancer cell HO-8910. The damping ratio of OVCAR-3 is close to that of HO-8910. The undamped frequency of ovarian cells decreases when the wave number increases, until it reaches a minimum value between 8 and 13, from which the undamped frequency starts increasing.

The undamped eigenfrequency of different continuum models is visualised in Figure 7. To plot this figure, the ovarian sample is assumed to be of the OVCAR-3 cell line. The nonlocal strain gradient theory (NSGT), traditional classical elasticity theory (TCET), strain gradient elasticity theory (SGET), and classical nonlocal elasticity theory (CNET) are taken into account. It can be concluded that the effect of the choice of the continuum model is more important at higher flexural wave numbers. The SGET predicts the highest undamped frequency while the eigenfrequency of the CNET is the lowest. The undamped eigenfrequency of the NSGT is lower than that of the TCET. In the NSGT, only the effect of strain gradient is considered, and both the first and second nonlocal parameters are set to zero. As nonlocal parameters are linked to structural softening, ignoring nonlocal effects leads to a significant increase in structural stiffness, leading to higher eigenfrequencies. Furthermore, strain gradients are strongly associated with stiffness hardening, thus further enhancing the tissue’s structural stiffness. The combined influence of stress nonlocality elimination and strain gradient results in a sharp increase in eigenfrequency in the SGET. Figure 8 shows the frequency parameter of flexural waves propagated within the ovarian samples when the dimensionless wave number is set to 20. The undamped frequencies of the ovarian cancer cells (both OVCAR-3 and HO-8910) are significantly lower than those of the healthy normal cells (HOSEpiC).

5.3. Breast Cancer

In the second case study, the frequency response of breast cancer cells to flexural waves is investigated in detail. Recent studies have demonstrated that the biomechanical properties of breast cells can be used as clues for breast cancer diagnosis and treatment efficacy assessments [50]. Two different biological conditions are considered: (1) MCF-7 human breast cancer cells; (2) MCF-10A healthy cells. The viscoelastic features of the breast cancer and healthy cells were obtained by Wang et al. [50] using scanning force microscopy. The average elasticity moduli of the MCF-7 and MCF-10A cells are 487.44 Pa and 1231.07 Pa, respectively. The structural viscoelasticity coefficient of the breast cancer cells (MCF-7) is 17.31 Pa·s, while the healthy breast cells (MCF-10A) have a viscoelasticity coefficient of 22.37 Pa·s. The Poisson’s ratio and mass density are set to 0.495 and 945 kg/m³, respectively. The biomechanical properties of the breast cells are assumed to be isotropic. The ratio of the side length of each breast sample to its thickness is set to 10, while the ratio of the sample length to its width is equal to one. The wave numbers in the x and y directions are assumed to be the same. Unless mentioned otherwise, a value of 0.05 is used for the dimensionless scale parameters. All breast cell samples are embedded in a viscoelastic medium made of dextran-based hydrogels with known viscoelastic properties. The elasticity modulus and viscoelastic coefficient of the surrounding medium are taken as 16 kPa and 4 kPa·s, respectively [54]. The ratio of the medium thickness to the sample thickness is 0.1.

In Figure 9, the imaginary part of the eigenfrequencies of MCF-10A and MCF-7 against the dimensionless flexural wave number is illustrated. The imaginary part of the eigenfrequency gradually decreases as the wave number increases until it reaches its critical value, where the imaginary part of the frequency vanishes. The healthy breast cells (MCF-10A) reach their critical value at the dimensional wave number 7.9. The critical wave number where the imaginary part of the eigenfrequency of breast cancer cells becomes zero is higher than that of the healthy cells at around 8.5. Figure 10 shows the real part of the eigenfrequency for the two breast samples versus the dimensionless wave number. It can be concluded that larger wave numbers are associated with greater differences between the real parts of the frequencies of the breast cancer and healthy cells. Before the critical point, there is only one curve for the real part of the eigenfrequency of the flexural waves, whereas there are two potential curves (branches) after the critical point. In both healthy and cancerous conditions, the lower branch is marked with an asterisk. The structural damping ratio and undamped frequency parameter are also visualised in Figure 11 and Figure 12 for the healthy breast cells (MCF-10A) and breast cancer cells (MCF-7), respectively. The damping ratios of both healthy and cancer cells enhance when the flexural wave number increases. Furthermore, the difference between the damping ratio of healthy and cancer cells is greater at higher wave numbers. For example, when the dimensionless wave number is set to 10, the damping ratio of MCF-10A is about 19% higher than that of the breast cancer cell MCF-7. Increasing the flexural wave number decreases the undamped frequency of breast cells until it reaches a minimum value between 11 and 16, depending on the pathological condition of the breast tissue. From the minimum point, the undamped frequency begins to increase with further increases in flexural wave number. The undamped eigenfrequencies of various continuum-based models are shown in Figure 13. The breast sample is assumed to be MCF-7. The influence of the underlying continuum model becomes more significant at higher flexural wave numbers. The SGET is linked to the highest undamped eigenfrequency, while the eigenfrequency of the CNET is the lowest. The NSGT estimates lower undamped eigenfrequencies than those predicted by the TCET. In Figure 14, the eigenfrequencies of the flexural waves propagated within the healthy breast cells and breast cancer cells are shown. The dimensionless flexural wave number is set to 20. The undamped eigenfrequency of the breast cancer cells (MCF-7) is lower than that of the healthy breast cells (MCF-10A).

5.4. Ovarian Fibrosis Associated with Ageing

The third case study is conducted on the flexural wave behaviour of ovarian fibrosis associated with ageing. Three different biological conditions are taken into consideration: (1) young mouse ovary, (2) old mouse ovary, and (3) old mouse ovary that has undergone collagenase treatment. The elasticity moduli of these ovaries are extracted from the available data reported in Ref. [51], which had been experimentally determined through an instrumental indentation technique. The mean elasticity constants of the young, old, and old collagenase ovaries are 1.98 kPa, 4.36 kPa, and 2.28 kPa, respectively. The structural viscoelasticity coefficient of the mouse ovaries is neglected in this case study. The Poisson’s ratio and mass density are set to 0.495 and 945 kg/m³, respectively. Unless stated otherwise, an amount of 0.05 is selected for all dimensionless scale parameters in the scale-dependent continuum model. The length-to-thickness ratio is set to 10 for all samples, whereas the aspect ratio of samples (i.e., length-to-width ratio) is equal to one. The wave number along the x axis is the same as the wave number in the y direction. The dimensionless definitions of the wave numbers in the x and y directions are given by $k_x{L}_x/\pi$ and $k_y{L}_y/\pi$, respectively. The effects of the Winkler, shear, and viscoelastic coefficients of the surrounding medium are not taken into account.

In the classical nonlocal model, the first nonlocal parameter is set to 0.05, while the second nonlocal parameter and the strain gradient parameter are set to zero. The eigenfrequency parameters of the mouse ovaries with different fibrosis conditions are shown in Figure 15a based on the classical nonlocal model. The frequency parameters of the young ovary are significantly lower than those of the old ovary. The difference is more pronounced at higher flexural wave numbers. It has been indicated that the collagenase treatment of the old mouse ovary has the potential to remarkably lower the eigenfrequency parameter, almost recovering the magnitude of the young ovary (Figure 16). Figure 15b illustrates the eigenfrequency of the different mouse ovaries, including old, young, and old with collagenase treatment, versus the flexural wave number for the strain gradient model. The nonlocal parameters are set to zero, while a value of 0.025 is used for the dimensionless strain gradient parameter. The frequency parameter of the young mouse ovary is the lowest, while the old ovary has the highest frequency parameter. The frequency of the old ovary that underwent collagenase treatment is close to that of the young mouse ovary, highlighting the importance of collagen in the frequency response of ovaries to transverse waves. Furthermore, the frequency responses of the old, young, and old collagenase ovaries under flexural waves for the traditional local model and the nonlocal strain gradient model are plotted in Figure 15c and Figure 15d, respectively. In the traditional local model, all scale parameters are neglected, while the nonlocal strain gradient model simultaneously incorporates the effects of the first nonlocal parameter (0.05), strain gradient (0.025), and the second nonlocal parameter (0.05). Compared to the classical nonlocal model, the frequency response predicted by the traditional local model is in a linear form, lacking the ability to mimic nonlinear responses.

5.5. Comparison Study of Different Machine Learning Models

In this sub-section, the results of various machine learning regressors are compared and discussed. These ML regressors include Ridge, Lasso, ElasticNet, SVM, and KNN, as well as random forest and the extreme gradient boosting technique. For more details about the ML approach and the training dataset, refer to Section 3 of this paper (machine learning approach).

Table 3. The optimised hyperparameters of different ML algorithms used for the prediction of the viscoelastic frequencies of biological tissues at diverse scales under flexural wave propagation.

Machine Learning Model	Best Hyper-Parameters	Best Cross- Validation Score	Test Mean Square Error
Random Forest	bootstrap: True, max_depth: None, min_samples_leaf: 1, min_samples_split: 2, n_estimators: 200	9.2463 × 10−6	4.2856 × 10−6
XGBoost	colsample_bytree: 0.7, gamma: 0, learning_rate: 0.1, max_depth: 3, min_child_weight: 3, n_estimators: 300, reg_alpha: 0.01, reg_lambda: 1.5, subsample: 0.8	4.9260 × 10−5	1.4257 × 10−5
Lasso	alpha: 0.1	0.00935	0.00829
Ridge	alpha: 0.1	0.00700	0.00642
ElasticNet	alpha: 0.1, l1_ratio: 0.1	0.00872	0.00770
SVR	C: 1, coef0: 0.5, degree: 4, epsilon: 0.01, gamma: scale, kernel: poly	4.8440 × 10−5	5.7263 × 10−5
KNN regressor	metric: Euclidean, n_neighbors: 3, weights: distance	9.1557 × 10−6	4.3665 × 10−6

shows the best parameters of each ML model obtained by the hyper-parameter tuning using the grid-search cross-validation technique. The ranking is conducted based on the negative mean squared error (MSE). The MSE of the test dataset is also listed in the table for the best estimator. It is found that linear regressors such as Lasso, Ridge, and ElasticNet fail to predict the frequency response of biological samples due to the complexity and intrinsic nonlinearity of the response. By contrast, the random forest and extreme gradient boosting models are capable of accurately predicting the frequency behaviour of biological samples under flexural wave propagation. In addition, the regressors of SVM and KNN give highly accurate estimations of the frequency response with a small MSE.

In

Table 4. The MSEs of the extreme gradient boosting algorithm developed for the flexural wave propagation analysis. Different numbers of estimators and various learning rates are considered.

XGBoost Learning Rate	Number of Estimators	Train Mean Square Error	Test Mean Square Error
0.01	5	0.00847	0.00757
	10	0.00776	0.00692
	100	0.00166	0.00143
	300	7.3643 × 10−5	6.5005 × 10−5
0.1	5	0.00376	0.00327
	10	0.00155	0.00131
	100	2.7657 × 10−6	1.4365 × 10−5
	300	2.5017 × 10−6	1.4257 × 10−5
0.2	5	0.00141	0.00122
	10	0.00024	0.00020
	100	5.5292 × 10−6	2.7804 × 10−5
	300	4.8781 × 10−6	2.6121 × 10−5

, the effects of two important parameters of the XGBoost model, namely the learning rate and number of estimators in the train and test MSEs, are investigated. It is found that the number of estimators plays a crucial role in the performance of the ML regressor. As the number of estimators increases, the MSEs of both the training and test data significantly decrease. Particularly, for the learning rate of 0.01, at least 300 estimators are required for a reasonable prediction, while 100 estimators are sufficient when the learning rate is either 0.1 or 0.2. The best estimation with the minimum MSE is achieved for a learning rate of 0.1 and 300 estimators. It should be noted that a higher number of estimators are associated with higher computational time, and thus, depending on the application and the required precision, a trade-off is considered between model accuracy and computational costs.

Table 5. Comparison between various ML regressors, including extreme gradient boosting (XGBoost), random forest, SVM, KNN, Ridge, and ElasticNet. The eigenfrequency is estimated for different randomly selected wave numbers between 0 and 20.

Wave Number	Actual Frequency	Machine Learning Model
		XGBoost	Random Forest	SVR	KNN	Ridge	ElasticNet
9.5	147.574	147.838	147.774	149.102	147.760	169.803	169.261
1.5	203.385	203.789	203.560	201.862	203.582	183.613	171.555
3	193.528	195.141	192.738	191.800	193.114	181.024	171.125
15.8	158.819	158.086	158.721	156.928	158.608	158.927	167.455
12.8	147.036	146.806	147.020	145.629	147.139	164.106	168.315
11.5	145.255	145.576	145.261	144.969	145.249	166.351	168.688
6.9	160.642	163.082	159.783	162.628	159.503	174.291	170.007
17.1	166.337	165.108	166.301	165.409	166.007	156.683	167.082
17.5	168.863	167.011	169.498	168.344	169.223	155.993	166.967
4.5	180.543	181.449	181.026	180.306	180.183	178.434	170.695
6.6	162.816	163.082	164.206	164.650	163.557	174.809	170.093
18.3	174.185	173.428	174.007	174.563	174.463	154.612	166.738
16.5	162.713	161.495	162.852	161.276	163.023	157.719	167.254
7.8	154.864	155.492	154.882	157.062	154.583	172.738	169.749
18.7	176.961	176.331	176.565	177.804	177.242	153.921	166.623

lists the eigenfrequencies of the ovarian cancer cell line OVCAR-3 for different wave numbers obtained by various machine learning techniques. The actual target values are also presented for comparison purposes. Fifteen different dimensionless wave numbers are randomly chosen between 0 and 20 to show the prediction capability of each model. It is observed that the results of the random forest technique and the KNN regressor closely match those of the actual data. A reasonable agreement is also found between the predictions made by the extreme gradient boosting algorithm and the actual data, and between the SVM regressor and the actual target frequencies.

6. Discussion

6.1. Ovarian Cancer

The imaginary part of the eigenfrequency is associated with damping, energy loss, and amplitude attenuation within the biological tissue during vibrations. The imaginary part can serve as an indicator of the viscoelastic damping coefficient. Any change in the patterns of the imaginary part of the eigenfrequency indicates alterations in the viscoelastic damping coefficient of the tissue. Recent experimental studies demonstrate that the tissue viscoelastic coefficient is significantly correlated with disease conditions, including cancer [55,56]. The visualisation of the imaginary part enables the detection of regions with abnormal damping properties, and consequently the identification of potential malignancy and diseased conditions. In this paper, both real and imaginary eigenfrequency parts have been plotted and discussed to allow for the comprehensive vibrational analysis of different biological tissues.

For all three different types of human ovarian cells, the imaginary part of the eigenfrequency decreases when the flexural wave number is increased (Figure 3). This is potentially associated with the fact that higher wave numbers lead to higher damping ratios, which gradually decrease the imaginary part, as can be seen from Equation (13). Healthy normal ovarian cells (HOSEpiC) reach their critical point, where their imaginary eigenfrequencies vanish, at lower flexural wave numbers than those of ovarian cancer cells (HO-8910 and OVCAR-3). This is associated with the faster growth of the damping ratio of healthy ovarian cells with an increase in the flexural wave number (Figure 5).

Another important finding is that the difference between the frequency responses of the ovarian cancer and healthy conditions is more pronounced at larger wave numbers (Figure 4). This enhanced difference indicates the significance of considering wave number in the detection of abnormalities based on biomechanical features. Areas of higher difference in flexural wave propagation behaviour would provide useful insights into the biomechanical response of ovarian cells, leading to more accurate detection methodologies. Strain gradient effects are associated with higher undamped eigenfrequencies, while stress nonlocality causes ovarian cells to behave in a softer manner (Figure 7). This is because strain gradients are linked to enhanced structural stiffness, while nonlocal effects lead to lower structural stiffness [20,28,30]. The flexural wave characteristics of the ovarian cancer cell line OVCAR-3 are similar to those of HO-8910 (Figure 6 and Figure 8), highlighting that these characteristics cannot be used to differentiate between different types of ovarian cancer. However, the undamped frequencies of ovarian cancer cells (both OVCAR-3 and HO-8910) are significantly lower than those of healthy normal ones (HOSEpiC) (refer to Figure 8), indicating that flexural wave characteristics could be used as diagnostic indicators at the cellular level for distinguishing between ovarian cancer and healthy ovarian tissue.

Overall, the proposed model has the potential to be used in biomechanical-based imaging techniques, such ultrasound elastography and MRI elastography, to detect tumours of different sizes. This capability is evidenced by incorporating the scale effects associated with cellular features. The incorporation of stress nonlocality and strain gradient allows for the accurate estimation of deformation behaviour at small-scale levels. However, these effects gradually decrease as the size of the tumour increases, and vanish after a certain size. For large tumours, the current model of viscoelasticity reduces to that of the classical model, which has been widely used for cancer detection [11,12]. Therefore, the present model could be helpful in the detection of solid tumours at different stages from early pre-malignancy to metastasis.

6.2. Breast Cancer

The imaginary part of the eigenfrequency decreases as the flexural wave number is enhanced for both breast cancer and normal healthy cells (Figure 9). This is because of the strong link between higher wave numbers and higher damping ratios. As mathematically demonstrated by Equation (13), the imaginary part of the eigenfrequency decreases with the increase in the structural damping ratio. Healthy normal breast cells (MCF-10A) reach their critical point of imaginary eigenfrequency at a lower flexural wave number than breast cancer cells (MCF-7). This is due to the faster growth of the damping ratio of healthy breast cells with increases in flexural wave number (Figure 11). Furthermore, it can be concluded that the difference between the frequency response of the breast cancer and healthy breast cells is enhanced when the flexural wave number increases (Figure 10 and Figure 11). This enhanced difference highlights the importance of wave number in elastography modelling for the detection of breast abnormalities. Regions of higher difference in flexural wave propagation behaviours would facilitate the detection of breast tumours at small-scale levels, such as microscales and cellular levels.

Strain gradient effects significantly enhance the undamped eigenfrequencies of transverse waves, while nonlocal stresses make clumps of breast cells behave in a softer manner (Figure 13). This is rooted in the fact that strain gradients are linked to enhanced structural stiffness, while nonlocal stresses yield lower structural stiffness [36]. The undamped frequencies of breast cancer cells (MCF-7) are considerably lower than those of healthy normal cells (MCF-10A) (Figure 12 and Figure 14). This finding further supports the fact that flexural wave characteristics at the cellular level can serve as an early indicator of breast cancer. The present non-classical viscoelasticity model contains three different scale parameters that enable the comprehensive incorporation of scale effects associated with cellular-level structures and features of biological tissues. This facilitates the development of scale-dependent elastography imaging techniques for the early detection of cancer by the visualisation of solid tumours at smaller scales, which are not visible with present imaging technologies. For more detail, refer to Table 1, which compares the capability of different continuum models used in elastography imaging to simulate scale effects.

6.3. Ovarian Fibrosis Associated with Ageing

In addition to ovarian and breast cancer detection, flexural wave characteristics can be used in the detection of ovarian fibrosis. The frequency parameters of the young mouse ovary are significantly lower than those of the old mouse ovary with increased fibrosis (Figure 15 and Figure 16). This is best evidenced by the fact that the old ovary has a higher amount of collagen, which consequently leads to higher structural stiffness. The difference is more pronounced at higher flexural wave numbers, highlighting the importance of the choice of wave number in the frequency behaviour of the ovary. Moreover, it has been shown that treating an aged mouse ovary with collagenase significantly lowers its eigenfrequency, restoring it to the magnitude observed in young mice. This highlights the importance of the collagen extracellular matrix in the frequency response of ovarian tissue subjected to flexural waves. Compared to the scale-dependent continuum models (NSGT, CNET, and SGT), the frequency response predicted by the traditional local model with no scale effects is in a linear form, indicating that the local continuum models are not able to describe complex nonlinear eigenfrequency responses (Figure 15). Nonlocal effects are associated with lower eigenfrequencies for all types of mouse ovaries, as stress nonlocality reduces structural stiffness [57]. By contrast, higher strain gradients lead to higher vibrational eigenfrequencies induced by the propagation of flexural waves, since there is strong correlation between strain gradient and structural stiffness at small-scale levels [20].

6.4. Comparison Study of Different Machine Learning Models

The application of linear regression models, including Lasso, Ridge, and ElasticNet, for the prediction of flexural wave propagation response in biomedical engineering and biomechanics is challenging (Table 3 and Table 5). The frequency responses of breast cancer, ovarian cancer, and ovarian tissue with enhanced fibrosis to flexural waves are inherently nonlinear at small-scale levels, making linear machine learning models ineffective. This finding is well aligned with the fact that biological systems often display nonlinear behaviours that are not fully captured through simple linear modelling techniques [58]. However, ML approaches such as random forest (Table 3 and Table 5) and extreme gradient boosting (Table 3, Table 4 and Table 5), with appropriate configuration and optimised hyper-parameters, are capable of the accurate prediction of the viscoelastic eigenfrequency responses of various biological samples, including healthy and cancerous ovary and breast lesions, as well as ovarian fibrosis. Furthermore, the SVM regressor with a polynomial kernel function and a gamma hyperparameter equal to ‘scale’ can make precise predictions with very small MSE values that are comparable to the random forest and extreme gradient boosting algorithms. The KNN regressor, with consideration of three neighbouring points and the weighting type of ‘distance’, also yielded highly accurate estimations of the viscoelastic eigenfrequency response (Table 3 and Table 5). However, the computational cost of this instance-based model could be high, as every time a new data point is analysed, the distances between the given points and all other points are required to be measured.

In ensemble algorithms such as XGBoost and random forest, the number of estimators plays a crucial role in the model’s precision and significantly affects the performance of the ML model [59]. Overall, it is concluded that higher numbers of estimators lead to lower MSEs in tasks such as the eigenfrequency analysis of human breast cancer, ovarian cancer, and even ovarian fibrosis by aging (Table 4). It should be noted that the choice of the number of estimators depends on the learning rate. For a learning rate of 0.01, a minimum of 300 estimators is needed to assure accurate forecasts. Nevertheless, 100 estimators are sufficient when higher learning rates of 0.1 and 0.2 are used. These findings demonstrate the importance of the hyperparameter adjustment process for obtaining a reasonable balance between the ML algorithms’ computational cost and precision. In scenarios where real-time predictions or rapid execution are required while maintaining computational resources, fewer number of estimators may be selected. By contrast, in applications focused on the accurate forecast of the flexural wave responses of breast and ovarian cells, a higher number of estimators is a sensible choice, despite the higher computational time and cost.

7. Conclusions

In this paper, a scale-dependent viscoelastic model has been developed to predict the flexural wave response of various pathological conditions, such as ovarian cancer, breast cancer, and ovarian fibrosis. The Kelvin–Voigt model of viscoelasticity, nonlocal strain gradient theory, and Hamilton’s law were used to derive the partial differential equations related to the time-dependent deformation of biological systems. The real and imaginary parts of the eigenfrequency, as well as the undamped frequency and the structural damping ratios, were obtained and analysed in detail. Different types of pathological conditions, including ovarian cancer, ovarian fibrosis associated with ageing, and breast cancer, as well as their healthy control counterparts, were taken into account. Several machine learning algorithms were developed and integrated with the nonlocal continuum approach, providing scalability, computational efficiency, and flexibility with the underlying assumptions of the continuum model. Using the grid-search cross-validation technique, hyper-parameter tuning was conducted to obtain the optimal configuration for each ML regressor. A comprehensive comparison study was performed on the SVM regressor and KNN regressor, on linear regression algorithms such as Lasso, Ridge, and ElasticNet, and on regressors of ensemble learning techniques like random forest and extreme gradient boosting.

The imaginary part of the eigenfrequency of both cancerous and healthy ovarian cells decreases when the flexural wave number increases, since higher wave numbers lead to higher damping ratios. HOSEpiC cells reach their critical point at lower flexural wave numbers than HO-8910 and OVCAR-3 cells. The difference between the eigenfrequency behaviour of ovarian cancer and healthy cells is more pronounced at larger wave numbers. Strain gradient effects are associated with higher undamped eigenfrequencies, while stress nonlocality causes ovarian cells to behave in a softer manner. Furthermore, the undamped frequencies of OVCAR-3 and HO-8910 cells are significantly lower than those of HOSEpiC cells, indicating that flexural wave characteristics can be used as diagnostic clues for the detection of ovarian cancer at the cellular level. In addition, it has been demonstrated that flexural wave characteristics can be utilised in the detection of other diseases such as breast cancer and ovarian fibrosis. Compared to the scale-dependent continuum models, the frequency response predicted by the traditional local model is almost linear, indicating that local models are not able to describe the complex eigenfrequency responses of biological systems. Furthermore, linear machine learning models fail to accurately predict the frequency responses of breast cancer, ovarian cancer, and ovarian fibrosis under flexural waves. However, machine learning approaches such as random forest and extreme gradient boosting with appropriate configuration and optimised hyper-parameters are highly effective in the flexural wave propagation analysis of biological systems.

8. Patents

Ali Farajpour, Wendy V. Ingman, “Scale-dependent elastography method for detection of small inclusions in biological tissue”. Applicant: The University of Adelaide, Application number: PCT/AU2023/050855. Patent Cooperation Treaty (PCT) Submission Date: 23 September 2023.