Pathological changes in the skin often vary the viscoelasticity of the skin^[1]. For example, systemic sclerosis makes the skin harder and thicker; the skin viscosity in patients with systemic sclerosis tends to grow with the degree of skin involvement and is approximately 1.5 to 3 times higher than that of healthy skin^[2]. Therefore, skin-viscoelasticity measurements can be a valuable tool in aiding the clinical diagnosis of skin diseases. Among the numerous viscoelasticity-measurement techniques, optical coherence elastography (OCE) is a noninvasive technique built on optical coherence tomography (OCT)^[3–5], which has the advantages of micron-level spatial resolution and sub-nanometer displacement sensitivity^[6]. OCE typically provides a quasi-static or dynamic measurement of elasticity; but quantifying the tissue viscoelasticity with OCE is in its infancy^[7,8]. The quasi-static OCE technique generally obtains the low-frequency viscoelasticity of tissues. The wave-based dynamic OCE technique can quantitatively measure tissue high-frequency viscoelasticity using elastic waves and is widely used in ophthalmology^[9], dermatology^[10], and other fields. In general, this technique uses a rheological model to fit the phase velocity of the elastic wave that propagates within a tissue under dynamic load excitation for obtaining the tissue viscoelasticity. Therefore, the choice of rheological models in this technique is extremely important for correctly calculating the viscoelasticity of tissues. A suitable rheological model should be able to accurately predict the mechanical properties of tissues over a wide range of frequencies and may provide specific parameters that serve as potential biomarkers in disease diagnosis.

In recent years, many rheological models have been proposed to correlate the elastic wave information with sample viscoelasticity. Rheological models with two to three parameters are commonly used in practical applications of elastography, which matches the general sense of Occam’s razor (where the simplest solution is preferred)^[11]. The Kelvin-Voigt (KV) model with two parameters is the most widely used model in OCE. Compared with the two-parameter Maxwell model, it is more consistent with the mechanical properties of soft biological tissues^[12]. However, the KV model is unable to predict the stress relaxation behavior, which is the time-dependent decrease of stress under constant strain. To solve this problem, researchers developed the three-parameter Kelvin-Voigt fractional derivative (KVFD) model by introducing fractional operators. More details on the application of fractional operators can be found in Ref. [13]. Compared with the KV model, the KVFD model can fit the stress relaxation curve of the tissue more accurately^[14,15] and reproduce the power-law behavior of cumulative multiple relaxation processes^[11]. The power-law behaviors occur frequently in tissue structures^[16].

Currently, many studies aim to obtain broadband elastic wave information for enhancing the reliability of viscoelastic calculations in OCE^[17,18]. Compared with low-frequency components of elastic waves, high-frequency components are more significant in determining the dispersive pattern or the viscous value^[12]. The absence of high-frequency information reduces the accuracy of viscoelasticity measurements in tissues such as the cornea^[19] and skin^[10]. However, to date, most studies on the performances of rheological models have been focused on low frequencies, and the performances of rheological models in calculating the high-frequency ($>1\mathrm{kHz}$) viscoelasticity of tissues have not been investigated adequately. Previous works have shown that the power-law behavior of living cells at high frequencies is different from that at low frequencies: the complex shear modulus follows a weak power law with an exponent between 0.05 and 0.35 at low frequencies (0.01–100 Hz); as the frequency increases, the complex shear modulus shows a stronger power law, and the loss modulus shows a greater dependence on frequency than the storage modulus^[20]. Thus, high-frequency storage and loss moduli of tissues affected by the high-frequency rheological behavior of living cells may have different frequency dependencies, which are rarely investigated in the performances of rheological models in OCE. The storage and loss modulus expressions of the KVFD model have the same power-law exponent, so it is not possible to accurately describe the storage and loss moduli with different dependencies on frequency.

To address this problem, a modified Kelvin-Voigt fractional derivative (mKVFD) model is constructed in this study based on the power-law behavior of soft tissues and the high-frequency complex shear modulus of living cells. The storage and loss moduli in the proposed model have independent dependencies on frequency. The performance of the mKVFD model was verified using the broadband (1 Hz–10 kHz) complex shear modulus of gelatin phantoms obtained through rheometer and OCE tests. Furthermore, the performance of the mKVFD model was also evaluated using the high-frequency (1–128 kHz) complex shear modulus of polydimethylsiloxane (PDMS) phantoms obtained from the open database. The test results demonstrate that the proposed mKVFD model can accurately represent the viscoelasticity of the measured samples over a wider range of frequencies.

The shear modulus of viscoelastic materials is a complex parameter $\widehat{G}(\omega )$ related to the complex Young’s modulus. The real and imaginary components of the complex shear modulus are the storage modulus $G_{\text{stor}}(\omega )$ and loss modulus $G_{\text{loss}}(\omega )$, respectively. The storage modulus is related to the elasticity of the tissue, and the loss modulus is related to the shear viscosity $\eta$ [where $\eta =G_{\text{loss}}(\omega )/\omega$]. When a shear wave propagates through a viscoelastic tissue, its complex wavenumber $\widehat{k}(\omega )$ exhibits the following relationship^[14]: $$\widehat{k}(\omega )=\omega \sqrt{\frac{\rho }{\widehat{G}(\omega )}}=\frac{\omega }{c_s(\omega )}-j\alpha (\omega ),$$(1)where $\rho$ is the tissue density, and $c_s(\omega )$ and $\alpha (\omega )$ are the phase velocity and attenuation coefficient of the shear wave, respectively. Therefore, the storage and loss moduli of tissues can be obtained from $c_s(\omega )$ and $\alpha (\omega )$ via $$G_{\text{stor}}(\omega )=\rho {\omega }^2{c}_s^2(\omega )\frac{[{\omega }^2-c_s^2(\omega ){\alpha }^2(\omega )]}{{[{\omega }^2+c_s^2(\omega ){\alpha }^2(\omega )]}^2},$$(2)$$G_{\text{loss}}(\omega )=\frac{2\rho {\omega }^3{c}_s^3(\omega )\alpha (\omega )}{{[{\omega }^2+c_s^2(\omega ){\alpha }^2(\omega )]}^2}.$$(3)The KV and KVFD models are commonly used in OCE to describe the viscoelasticity of tissue-mimicking materials. In the KV model, the storage modulus of the sample is estimated using a constant shear modulus ${\mu }_0$, and the loss modulus is estimated using the equation $G_{\text{loss}}=\omega {\eta }_0$, where $\omega$ is the angular frequency and ${\eta }_0$ is the shear viscosity^[21]. The storage and loss moduli of the KVFD model in the frequency domain are, respectively, defined by the following equations^[15]: $$G_{\text{stor}}(\omega )={\mu }_0+{\mu }_a{\omega }^a\cos \left(\frac{\pi }{2}a\right),$$(4)$$G_{\text{loss}}(\omega )={\mu }_a{\omega }^a\sin \left(\frac{\pi }{2}a\right),$$(5)where ${\mu }_a$ is the viscous parameter and $a$ is the power-law exponent. The power-law exponent with a value between (0,1) is physically interpreted as the fluidity of tissues^[22]. In particular, the KVFD model can be simplified to the KV model when $a=1$. From Eqs. (4) and (5), the two power-law functions have the same power-law exponent. Thus, the storage and loss moduli of the KVFD model have the same dependence on frequency.

Considering the power-law rheological behavior in living cells at high frequencies (i.e., the loss modulus shows a greater dependence on frequency than the storage modulus)^[20,23], a modified KVFD model is constructed in this study. The frequency-domain expression of the modified rheological model is given by $${\widehat{G}}_{\mathrm{mKVFD}}(\omega )={\mu }_0+{\mu }_1{\omega }^{a_1}+j{\mu }_2{\omega }^{a_2},$$(6)where ${\mu }_1$ is the elastic parameter, $a_1$ is the elastic power-law exponent, ${\mu }_2$ is the viscous parameter, and $a_2$ is the viscous power-law exponent. The first term of this expression describes the elastic response of permanent structures in tissues such as cell walls; the second and third terms describe the power-law behavior of tissue structures such as the cytoskeleton. Power-law behavior is common in biological multiscale systems. The two power-law functions with different power-law exponents describe the storage and loss moduli in this model, respectively. Thus, the storage and loss moduli of this model have independent dependence on frequency.

In particular, the mKVFD model can be simplified to the KVFD model when $a_1=a_2$; the mKVFD model can be simplified to the KV model when $a_1=0$ and $a_2=1$. For viscoelastic samples with a negligible shear modulus ${\mu }_0$, the mKVFD model reduces to the dual power-law model [${\widehat{G}}_{\mathrm{dp}}(\omega )={\mu }_1{\omega }^{a_1}+j{\mu }_2{\omega }^{a_2}$].

Recent studies have shown that compared with the KV model, the KVFD model can be simplified to a spring-pot model [${\widehat{G}}_{\mathrm{SP}}(\omega )={\mu }_a{(j\omega )}^a$], which can more accurately characterize the stress relaxation curve of liver tissue in the time domain and the elastic wave dispersion curve in the frequency domain (1 Hz–1 kHz)^[24]. However, the loss factor (i.e., the ratio of loss to storage modulus) in the spring-pot model is constant, whose value is equal to $\tan (\pi a/2)$, and this does not match the frequency-dependent increase of the loss factor measured in tissues^[25,26]. On the contrary, the loss factor in the dual power-law model is not constant and is related to ${\omega }^{a_2-a_1}$. Therefore, compared with the spring-pot model, the dual power-law model can better describe the loss factor in biological tissues.

In this study, gelatin (a commonly used tissue-mimicking material) phantoms were used as experimental samples to evaluate the performance of the proposed model and other rheological models. First, to evaluate the performance of the mKVFD model in the low-frequency range, the storage and loss moduli of gelatin phantoms in the frequency range of 1–50 Hz were measured using a commercial rheometer. Subsequently, the high-frequency complex shear moduli of gelatin phantoms were obtained using the OCE system under pure tone excitation. Because of the crucial importance of high-frequency surface wave information at 4–10 kHz for skin-elasticity calculations^[10], the maximum excitation frequency of surface waves in the OCE test was set at 10 kHz. The minimum excitation frequency was selected as 1 kHz. The complex shear modulus of PDMS phantoms with a frequency range of 1–128 kHz from the open database was added to further validate the accuracy of the mKVFD model in high-frequency viscoelasticity calculations of tissue phantoms other than gelatin.

In this study, pig skin gelatin (Type A, 300Bloom, Sigma-Aldrich Corp., USA) phantoms with different mass concentrations (10%, 12%, 14%, and 16%) were used to simulate tissues with different viscoelasticities. Meanwhile, 0.6% (mass concentration) titanium dioxide (20 nm particle size, Marklin, China) was added to the gelatin solution. This was performed to increase the scattering of samples. In this study, the rheometer test used samples with a diameter of 25 mm and a thickness of approximately 2 mm; additionally, the size of samples used for the OCE test was approximately $40\mathrm{mm}(\mathrm{diameter})\times 12\mathrm{mm}(\mathrm{height})$. The elasticity of gelatin phantoms is susceptible to temperature^[27]; thus, all experiments in this study were performed at the constant room temperature (22.5°C). To measure the phantom densities, the mass of each phantom was obtained using a standard laboratory scale (JNB6002, JOANLAB, China), and the volume of each phantom was measured using the water displacement method^[28]; finally, the average density of the gelatin phantoms used in the experiments was calculated to be $1038\mathrm{kg}/{\mathrm{m}}^3$.

Rheometer tests were performed using a commercial shear rheometer (KINEXUS lab+, NETZSCH, Germany) with a 25 mm diameter parallel plate configuration. The sample temperature was stabilized to the temperature (22.5°C) set by the temperature control unit before the commencement of the test. The linear viscoelastic region (LVR) is the strain range in which the complex shear modulus remains constant. First, strain amplitude-sweep tests with the strain amplitude incrementally increasing from 0.01% to 10% were performed on the gelatin phantoms at a frequency of 1 Hz to determine the LVR. Based on this, a shear stress of 25 Pa (corresponding to less than or equal to 0.5% of the complex shear strain) was selected within the LVR to perform subsequent frequency-sweep tests in the frequency range of 1–50 Hz. The storage and loss moduli of gelatin phantoms with four different concentrations shown in Sec. 3.1 (two samples for each concentration) were then obtained.

The OCE system used in this study comprised an elastic wave detection subsystem and an excitation subsystem, as depicted in Fig. 1. The elastic wave detection subsystem was a swept-source OCT system that tracked the propagation of elastic waves in the M-B scan mode. Details of this swept-source OCT system and the M-B scan mode were presented in the previous works^[29]. Briefly, this system included a laser (SANTEC, Japan) with a center wavelength of 1300 nm, a bandwidth of 100 nm, and an A-scan rate of 50 kHz. The lateral resolution of the system was 14.5 µm, and the axial resolution of the system was 11.2 µm. The axial displacement $u_z(z,x,t)$ of the given position $x$ can be obtained using the following equation: $$u_z(z,x,t)={\int }_{t_1}^{t_2}\frac{u_{z,\tau }(z,x,t)}{\tau }\mathrm{d}t={\int }_{t_1}^{t_2}\frac{{\lambda }_0\mathrm{\Delta }\phi (z,x,t)}{4\pi n\tau }\mathrm{d}t,$$(7)where $\tau$ is the adjacent A-line time interval, $u_{z,\tau }(z,x,t)$ is the axial displacement at a given position between two adjacent A-lines, ${\lambda }_0$ is the central wavelength of the light source, $n$ is the refractive index of the air above the sample when measuring the displacement of the sample surface, and $\mathrm{\Delta }\phi (z,x,t)$ is the Doppler phase shift between two adjacent A-lines. Trigger signals of the laser were used to control the synchronization of the data acquisition and the excitation trigger signal. At each transverse location, an M-scan acquired for each excitation was composed of 512 consecutive A-lines, and 512 sampling positions along the lateral scanning direction were recorded to form a complete M-B scan.

In the elastic wave excitation subsystem, a customized contact probe bonded to a piezoelectric transducer was used for the excitation. Considering that the linear wave has lower geometrical attenuation than the circular wave^[10], this probe was designed to have a 2 mm contact edge with the sample to produce elastic waves with linear wavefronts. When generating a peak-to-peak signal, a pure tone (sine) excitation can produce the maximum power at an individual frequency compared with pulsed excitation. Therefore, a signal generator (DG4062, RIGOL, China) was used to generate a pure tone excitation waveform in the external trigger mode, which was amplified by a power amplifier to excite surface waves on the sample surface. In this case, the frequency range of the pure tone excitation mode was 1–10 kHz with an interval of 1 kHz.

The coefficient of determination (also known as $R^2$) is a statistical measure that indicates the degree of fit between the data points and the regression curve^[30]. An $R^2$ value closer to 1 indicates the higher closeness of the data points and the fitted curve. In this study, $R^2$ was used to evaluate the goodness-of-fit between the model and the test values, which can be calculated using the following equation: $$R^2=1-\frac{\sum _{i=1}^n{(y_i-f_i)}^2}{\sum _{i=1}^n{(y_i-\overline{y})}^2},$$(8)where $\overline{y}$ is the mean value of the test data $y_i$ and $f_i$ is the value on the regression line. However, the $f_i$ of the storage modulus in the KV model equals the $\overline{y}$; thus, the goodness-of-fit concerning the storage modulus cannot be evaluated using $R^2$. The root mean square error ($\mathrm{RMSE}=\sqrt{[\sum _{i=1}^n{(y_i-f_i)}^2]/n}$) was used to evaluate the degree of fit between the storage modulus of samples and the model parameters in this study.

The storage and loss moduli of the gelatin phantoms with four different concentrations were obtained from the rheometer test, as depicted in Fig. 2. The error bars in Fig. 2 represent the standard deviation at each frequency. In the rheometer test, the high-frequency data (about $>20\mathrm{Hz}$) were strongly affected by the inertial effects of the sample and the rheometer^[31] and were prompted by rheometer equipment as unreliable data; thus, these high-frequency data were not used in the curve fitting and were not plotted in Fig. 2. Figure 2 shows that the storage modulus varies slightly in the frequency range of 1–20 Hz, and the loss modulus increases with increasing frequency.

By substituting the storage and loss moduli into the KV model, the viscoelastic parameters of the KV model (i.e., ${\mu }_0$ and ${\eta }_0$) can be obtained. The ${\mu }_0$ was estimated by averaging the storage modulus at different frequencies; the linear regression of the loss modulus was calculated to obtain the ${\eta }_0$^[21,32]. The curves depicted in Figs. 2(a) and 2(d) show the fitting results of the storage and loss moduli using the KV model. The fitting parameters of the KV model (mean ± standard deviation) are listed in Table 1. Both ${\mu }_0$ and ${\eta }_0$ increase with the gelatin concentration, which is consistent with previous reports^[33,34]. In addition, the $R^2$ values of fitting the loss modulus using the KV model in all gelatin phantoms are less than 0, which indicates the failure of the KV model in predicting the variation of the loss modulus with frequency.

For the viscoelastic parameters of the KVFD model, the loss modulus expression [Eq. (5)] was fitted to the loss modulus data obtained by the rheometer using the least-squares nonlinear regression method to estimate the ${\mu }_a$ and $a$. The obtained ${\mu }_a$ and $a$ values were substituted into the storage modulus expression [Eq. (4)] as known quantities, which were used to fit the storage modulus data, thereby estimating the ${\mu }_0$. The results of fitting the rheometer data using the KVFD model are shown in Figs. 2(b) and 2(e). The 95% confidence interval represented by the shaded area in the KVFD model was smaller than that present in the KV model. The $R^2$ values of fitting the loss modulus using the KVFD model range from 0.85 to 0.95, as listed in Table 1. The $R^2$ value of the KVFD model was higher than that of the KV model for the gelatin phantom with each concentration. However, for the storage modulus, the RMSE values of the KVFD model were slightly higher than those of the KV model.

The viscoelastic parameters of the proposed mKVFD model were estimated using the least-squares nonlinear regression method. The curves shown in Figs. 2(c) and 2(f) demonstrate the fitting results of the mKVFD model for the storage and loss moduli measured using the rheometer, respectively, and the fitting parameters are shown in Table 1. The mKVFD model obtains the same fitting results as the KV model for the storage modulus and the same goodness-of-fit as the KVFD model for the loss modulus. Therefore, the mKVFD model is superior to the KVFD model in the characterization of the low-frequency storage modulus and has the same advantages as the KVFD model in the characterization of the low-frequency loss modulus.

In addition to the storage and loss moduli of the samples, the rheometer was also able to obtain shear viscosity $\eta$. The $\eta$ of gelatin phantoms at different concentrations were fitted using linear regression in logarithmic coordinates with $R^2$ values greater than 0.99, as shown in Fig. 3; thus, $\lg \eta$ and $\lg \omega$ have a linear relationship, which is a power-law relationship with $\eta$ and $\omega$. This power-law relationship can be derived from the mKVFD model [$\eta =G_{\text{loss}}(\omega )/\omega ={\mu }_2{\omega }^{a_2-1}$], which further validates the accuracy of the loss modulus expression in the improved model.

To verify the performance of the mKVFD model in high-frequency viscoelastic calculations, the M-B scanning data of the gelatin phantoms were obtained using the OCE system to calculate the phase-velocity dispersion curves of the high-frequency surface wave. First, the intensity gradient was utilized to identify the sample surface, and the axial displacement $u_{z,\tau }(z,x,t)$ of the sample surface was obtained according to Eq. (7). Then, a two-dimensional Fourier transform was performed on the spatial-temporal slice of the displacement field to obtain the wavenumber-frequency domain. Finally, the maximum wavenumber of the surface wave at the excitation frequency was determined to obtain the phase velocity of the surface wave $c_R(f)=2\pi f/k$. The attenuation of the surface wave amplitude increases with increasing frequency. The maximum frequency of the phase-velocity dispersion curve was determined by the criterion that one full wave should be observed in the sample. According to the obtained phase velocity, the maximum wavelength of surface waves (2.00–2.93 mm), which is approximately equal to the maximum penetration depth of surface waves, is less than the thickness of the sample (12 mm). It relatively satisfies the infinite depth assumption for the Rayleigh surface wave. In this case, the frequency dependence of surface wave velocities is attributed to the material dispersion^[19].

The surface-wave attenuation coefficients of gelatin phantoms were obtained from the wavenumber-frequency domain using the AMUSE method^[35] to obtain the high-frequency complex shear modulus of gelatin phantoms. The surface-wave attenuation coefficients were calculated by measuring the full-width at half-maximum (FWHM) of each frequency peak $\mathrm{\Delta }{k}_{\mathrm{FWHM}}$ in the wavenumber curves using the following equation: $$\alpha =\frac{\mathrm{\Delta }{k}_{\mathrm{FWHM}}\times \pi }{\sqrt{3}}.$$(9)

Based on the surface wave velocities and attenuation coefficients, the storage and loss moduli of the gelatin phantoms with four concentrations were calculated and then fitted using the KV, KVFD, and mKVFD models, respectively, to obtain the results shown in Fig. 4.

In the KV model, a constant shear modulus was used to fit the high-frequency storage modulus that increases with frequency, leading to significant errors. Compared with the fitting results of the KV and KVFD models, the fitted curves of the mKVFD model are closer to the measured values of the storage modulus, especially in the gelatin phantoms with high concentration. The 95% confidence interval represented by the shaded area in the mKVFD model was significantly smaller than that present in the KVFD model. Table 2 lists the fitting parameters of the KV, KVFD, and KVFD models. The fitting error of the KV model for the storage modulus ($\mathrm{RMSE}=1.04–4.91$) is approximately twice that of the KVFD model ($\mathrm{RMSE}=0.391–2.70$, $R^2=0.70–0.86$). The mKVFD model provides an excellent fit to a high-frequency complex shear modulus with the value of ${\mu }_0$ being negligible, resulting in a dual power-law model with accurate results. The power-law exponents of the mKVFD model in terms of the storage modulus are different from those in terms of the loss modulus, where the model fitting errors ($\mathrm{RMSE}=0.134–0.843$, $R^2=0.96–0.98$) are further reduced. Thus, compared with the KV and KVFD models, the storage modulus expression of the mKVFD model is more accurate. The results of fitting the loss modulus with the mKVFD model are consistent with that with the KVFD model, and the goodness-of-fit of the mKVFD model is higher than that of the KV model. Overall, the mKVFD model is more appropriate for describing the high-frequency mechanical behavior of gelatin phantoms.

The viscoelastic properties of gelatin phantoms measured from the rheometer and OCE tests may vary due to factors associated with the test conditions, such as the stress magnitude and preloading effect^[24]. Therefore, the low-frequency complex shear modulus obtained from the rheometer test and the high-frequency complex shear modulus obtained from the OCE test were separated for model fitting and analysis in this study. The experimental results show that the power-law exponents of the modified model obtained by fitting the high-frequency complex shear modulus of gelatin phantoms are higher than those obtained by fitting the low-frequency complex shear modulus. Another possible reason is that gelatin phantoms, like cells, are viscoelastic materials exhibiting diverse mechanical behaviors at different timescales. The complex shear modulus of gelatin phantoms at high frequencies has a greater dependence on frequency than that at low frequencies.

The PDMS is also a skin-mimicking material that is commonly used in OCE^[10,36]. This study further uses the mKVFD model to fit the high-frequency (1–128 kHz) complex shear modulus of the PDMS from publicly available data in Ref. [18]. Figure 5 shows the results of the model fitting the high-frequency complex shear modulus of the PDMS. The high-frequency storage and loss moduli of the PDMS increase with frequency. For the storage modulus, the shear modulus of zero frequency ${\mu }_0$ was 0.65 MPa, the elastic parameter ${\mu }_1$ was $80.47\mathrm{Pa}·{\mathrm{s}}^{{\mathrm{a}}_1}$, and the elastic power-law exponent ${\mathrm{a}}_1$ was 0.79 in the mKVFD model. The fitting error of the mKVFD model ($\mathrm{RMSE}=0.19\mathrm{MPa}$, $R^2=0.98$) is approximately 83% lower than that of the KV model ($\mathrm{RMSE}=1.15\mathrm{MPa}$), which describes the storage modulus using a constant. For the loss modulus, the viscous parameter ${\mu }_2$ was $3.32\mathrm{Pa}·{\mathrm{s}}^{{\mathrm{a}}_2}$, and the viscous power-law exponent ${\mathrm{a}}_2$ was 1.00, with a goodness-of-fit greater than 0.99 in the mKVFD model. It is evident that the loss modulus of the PDMS has a higher frequency dependence than the storage modulus at high frequencies. This is similar to the cellular high-frequency power law behavior. For the KVFD model, the power-law exponent of 1.00 is obtained when fitting the loss modulus, at which point the KVFD model simplifies to the KV model, resulting in the model obtaining the same fit result as KV for the storage modulus. On the contrary, the mKVFD modulus with independent power-law exponents on the storage and loss moduli can accurately fit the high-frequency storage modulus. Therefore, the mKVFD model proposed in this study possesses the potential to predict the mechanical behavior of a wide range of tissue-mimicking materials and soft tissues.

A modified KVFD model based on the power-law behavior of soft tissues and the dependence of the high-frequency shear modulus in living cells on frequency is proposed in this study, and it has different power-law exponents for the storage and loss moduli. Thus, the storage and loss moduli in this modified model have independent dependencies on frequency, which can solve the problem of the spring-pot model not being able to accurately predict the loss factor variation with frequency in soft tissues. For the low-frequency complex shear modulus obtained from the rheometer test, the modified KVFD model can not only describe the approximately constant low-frequency storage modulus of gelatin phantoms but also provide the prediction of the power-law relationship between the low-frequency shear viscosity and frequency ($\mathrm{accuracy}>99\%$). To obtain the high-frequency shear modulus of gelatin phantoms, the surface wave dispersion curves and attenuation coefficients in the frequency range of 1–10 kHz were obtained by the constructed OCE system in this study. The results of the OCE test show that, compared with the KV and KVFD models, the modified KVFD model has a higher goodness-of-fit ($\mathrm{accuracy}>96\%$) for the high-frequency storage modulus of gelatin phantoms. To further validate the performance of the modified model in calculating the viscoelasticity of multiple tissue-mimicking phantoms, the model performance was tested using the complex shear modulus of the PDMS in the frequency range of 1–128 kHz from publicly available data. The results show that the modified KVFD model is also able to accurately fit the high-frequency storage modulus ($R^2=0.98$) and loss modulus ($R^2=0.99$) of the PDMS. In conclusion, the proposed modified model has the potential to accurately predict the mechanical properties of soft tissues over a wider range of frequencies, thereby providing insights into structural changes within biological tissues by monitoring the corresponding parameters of that model.