# A Mathematical Model to Study the Mechanical Information Induced by Lifting-Thrusting Needle

Authors: Yi Yu; Wei Yao; Guanghong Ding
Journal: Evidence-Based Complementary and Alternative Medicine
 (2019)

---

## Abstract

Focusing on the mechanical effect of traditional Chinese acupuncture, this study builds a mathematical model that simulates the mechanical process of lifting-thrusting needle. Analytic and numerical solutions are obtained to explore the mechanical information (displacement, strain, stress, and energy) in the skin tissue. Our results show that(1) needle manipulation leads to tissue displacement and mechanical stress field, but the needle should be inserted into the right position (about π/ω cm around the acupoint, where ω is the angular frequency) and enough depth (about 2 cm in lower limbs) to achieve effective mechanical stimuli; (2) the tissue displacement decays with an increase of distance from the stimulus position, more rapidly at higher frequencies; (3) there is an inverse relationship between the area of the ‘effective influence region’ (where shear strain > 0.2) and the stimulus frequency, which means larger needle movement is needed at higher frequencies to achieve a better curative effect; (4) more energy is required to maintain high frequency manipulation. This study proposes a preliminary comprehension of the mechanical response around the needle during the acupuncture process.


---

## Body

## 1. Introduction

Acupuncture, a physiotherapy sourced from Traditional Chinese Medicine (TCM), now gains recognition all around the world. Its therapeutic effect has been widely accepted since the National Institutes of Health Hearing in 1997 [1]. Acupuncture is a biomechanical process, during which a needle is inserted into the skin at the acupoint, and mechanical stimulation to the tissue is generated when twisting or lifting the needle. Elucidating the mechanical and biological mechanism of acupuncture effect is one of the urgent problems in TCM as well as international researches on physiotherapy.The premise to accomplish acupuncture effect is to get sensations (“De Qi” in Chinese), which means to maintain a certain intensity of stimulation by manipulating the needle. However, “De Qi” is usually described subjectively according to the experiences of skillful acupuncturists and the feelings of some sensitive patients; the underlying issues are lack of quantitative research. Therefore, it is important to make clear the mechanical sensation in acupuncture and to reveal the key factors related to it.Numerous work and different explanations have been put forward to help understand the acupuncture process. Recent morphological studies have made clear that the acupoint is mainly composed of connective tissue. As shown in Figure1, the dermic dense connective tissue and subcutaneous loose connective tissue form a three-dimensional collagen fiber network with numerous intertwining blood vessels, nerves, mast cells, and lymphatic vessels [2, 3]. Goh et al. hypothesized that there might be a relationship between the efficacy of acupuncture treatment and the depth of insertion needle [4]. Langevin et al. hypothesized that needle manipulation generates mechanical signals at acupoints [5], where longitudinal waves and transverse waves can propagate. Acupuncture needle manipulation results in sustained stretching by lifting-thrusting or rotating in the tissue of the acupoint [6]. Yu et al.’s experiment showed that the acupuncture effect and the force on the needle are very significant when the collagen structure is complete, and when the collagen is destroyed, the acupuncture effect is significantly weakened and the force on the needle body is very small [7]. Mvogo et al. constructed the analytic solutions of the tissue with collagen [8].Figure 1
Acupoint morphology. (a) Basic pattern of organization of subcutaneous layers [9, 10]; (b) acupuncture location in lower limbs.
(a) (b)The mechanical effect of needle insertion and manipulation is the key to acupuncture; however, researches related with the quantitative mechanical effect of acupuncture are limited. Deleuze et al. has numerically simulated the features of subcutaneous interstitial flow induced by acupuncture [11]. Thiriet et al.’s mathematical model analyzed the different effects between acupuncture positioned at and out of acupoints [12]. We have also constructed a mathematical model to analyze the needle rotation at acupoints [13]. In this study, the previous work is extended, and the acupuncture of lifting-thrusting manipulation is studied through the establishment of a mechanical model based on the theory of viscoelastic mechanics. The influences of needle factors, such as stimulus frequency, movement amplitude, radius, and insertion depth, have been discussed in order to explore the acupuncture mechanical effect.


## 2. Methods

### 2.1. Model

Soft tissues are usually assumed as viscoelastic materials; some are more “elastic” and some are more “viscous”; they exhibit obvious hysteresis, stress relaxation, and other viscoelastic phenomena [14]. In our model, the skin tissue is treated as a semi-infinite viscoelastic body and the needle is treated as a rigid body. As shown in Figure 2, the tissue is divided into two parts. Part (1) represents the tissue from the plane of the tip location to infinity (z ≥0) where tissue displacement is caused by the movement of the needle tip (point A). Part (2) is the needle manipulation region (– L0≤ z < 0), where L0 is the needle insertion depth and the standard L0 is 2 cm in lower limbs.Figure 2
Diagrammatic sketch of lifting-thrusting manipulation.


### 2.2. Equations

The soft tissue is supposed as a homogeneous isotropic incompressible viscoelastic body and fit for Feng's quasilinear viscoelastic theory of small amplitude displacement [15]. For acupuncture manipulation, the needle movement causes small amplitude oscillatory displacements in the low stress range (physiological state) of the tissue. The boundary conditions are specified as simple harmonic functions of time; therefore, the variables can be expressed as(1)uixi,t=u~ieiωtσijxi,t=σ~ijeiωtεijxi,t=ε~ijeiωti,j=1,2,3.where ui is the displacement component which is the distance from the original (equilibrium) position, xi is the coordinate component, t is the time, ω is the angular frequency, σij is the component of stress tensor, and εij is the component of strain tensor.Based on the strain-displacement relation and excluding the effect of volume force, we get the governing equation of motion in linear viscoelastic medium(2)λ∗iω+μ∗iωduj,ji+μ∗iωdui,jj=ρ∂2ui∂t2.where μ∗(iω) and λ∗(iω) are the viscoelastic material parameters andρ is the density. Postulate the tissue is incompressible, that is, uj,ji=0, and putting (1) into (2), we get(3)μ∗iωdu~i,jj=-ρω2u~iRewrite it in tensor form(4)μ∗iω∇2u~→=-ρω2u~→Suppose the tissue around the needle is a semi-infinite space. This is an axisymmetric problem (uθ=0,∂/∂θ=0). Replace μ∗ with μ and u~ with u; then we have (5) in cylindrical coordinate(5)μ∂2ur∂r2+1r∂ur∂r-urr2+∂2ur∂z2=-ρω2urμ∂2uz∂r2+1r∂uz∂r+∂2uz∂z2=-ρω2uz.where ur and uz are the displacement amplitude inr andz direction, respectively. Defining the dimensionless parameters: UZ=uz/u0,UR=ur/u0,Z=z/u0,R=r/u0, where u0 is the amplitude of needle movement. Then (5) is changed to(6a)μ∂2UR∂R2+1R∂UR∂R-URR2+∂2UR∂Z2=-ρω2u02UR(6b)μ∂2UZ∂R2+1R∂UZ∂R+∂2UZ∂Z2=-ρω2u02UZ.Based on the strain-displacement relation(7)εrz=12∂uz∂r+∂ur∂zeiωt=12∂UZ∂R+∂UR∂Zeiωt=ERZeiωt.where ERZ is calculated by explicit difference scheme(8)ERZ=12∂UZ∂R+∂UR∂Z≈12ΔRUZj+1,k-UZj,k+12ΔZURj,k+1-URj,k.And the corresponding stress is(9)τrz=Re2μiωεrz=Reμ1+iμ2·2ERZcos⁡ωt+isinωt=2ERZμ1cos⁡ωt-μ2sin⁡ωt=2ERZμcos⁡ωt+δ.where μ=μ12+μ22, tan⁡δ=μ2/μ1, and δ is the phase difference between stress and strain. The shear stress amplitude TRZ=2μERZ.In a cycle, the energy dissipation of one-unit volume caused by lifting-thrusting is(10)ΔE=∫02π/ωτrzε˙rzdt=4μERZ2∫02π/ω-ωμ1cos⁡ωt-μ2sin⁡ωtsin⁡ωtdt=4πERZ2μ2.The energy dissipation of the whole domain is(11)Ef=∫r0/u0+∞∫-L0/u0+∞ΔEdZdR=∫r0/u0+∞∫-L0/u0+∞4πμ2ERZ2dZdR.Define Emax as the max energy stored(12)Emax=4πERZ2∫0π/ωμ1ωcos⁡ωtsin⁡ωtdt=2πμ1ERZ2.Therefore, the mechanical quality factorQ which reflect the ratio of energy storage to energy dissipation is(13)Q=2πEmaxΔE=μ1μ2.


### 2.3. Boundary Conditions

The motion of the boundary is(14)ur=0,uz=u0cos⁡ωtr=r0,-L0≤z≤0ur=uz=0z→∞uz=0r→∞,wherer0 is the radius of the needle. Then(15)UR=0R→0∪R→∞∪Z→∞UZ=0R→∞∪Z→∞UZ=1R→0∩-L0u0≤Z≤0.


### 2.4. Analytic Solutions of Part (1)

LetUR=F(R)G(Z); (6a) is expressed as(16)FRR+1/RFR-1/R2FF+GZZG+ρω2u02μ=0.Let GZZ/G=λ; then we get(17a)GZZ-λG=0(17b)R2FRR+RFR+ρω2u02μ+λR2-1F=0.Equation (17b) is a 1-order Bessel equation, and the solutions are(18)GZ=A1eλZ+A2e-λZFR=B1J1ρω2u02μ+λR+B2Y1ρω2u02μ+λR,whereJ1,Y1 are 1-order Bessel function and Neumann function, respectively.InputtingUR(Z→∞)=0, we get A1=0. Because J1(0)=0 and Y1(0)=-∞, based on the boundary condition UR(R→0,Z=0)=0, we get B2=0; therefore(19)UR=AJ1ρω2u02μ+λRe-λZ.Because J1(+∞)=0, then UR(R→∞)=0.PutUZ=F(R)G(Z) into (6b), (20)FRR+1/RFRF+ρω2u02μ+GZZG=0.Let GZZ/G=λ; then we get(21a)GZZ-λG=0(21b)FRR+1RFR+ρω2u02μ+λF=0.Equation (21b) is a 0-order Bessel equation; therefore the solutions are(22)GZ=A1eλZ+A2e-λZFR=B1J0ρω2u02μ+λR+B2Y0ρω2u02μ+λR.InputtingUz(Z→∞)=0, we get A1=0. Because J0(0)=1 and Y1(0)=-∞, based on the boundary condition Uz(R→0,Z=0)=1, we get B2=0; therefore(23)UZ=J0ρω2u02μ+λRe-λZ.Suppose the tissue is incompressible, that is,∇∙v→=0; therefore(24)∂UR∂R+URR+∂UZ∂Z=0.Putting (19) and (23) into (24), we get(25)A∂J1ρω2u02/μ+λR∂R+AJ1ρω2u02/μ+λRR-λJ0ρω2u02μ+λR=0.Based on the recurrence relation of Bessel’s function(26)Jv-1x+Jv+1x=2vxJvxJv-1x-Jv+1x=2Jv′x,we get (27)∂J1kR∂R=kJ0kR-J2kR2J1kRR=kJ0kR+J2kR2k=ρω2u02μ+λ.Putting (27) into (25), we get (28)AkJ0kR-J2kR2+J0kR+J2kR2-λJ0kR=0,and then A=λ/k. Therefore(29)UZ=J0ρω2u02μ+λRe-λZUR=λρω2u02/μ+λJ1ρω2u02μ+λRe-λZ.


### 2.5. Calculation Solutions of Part (2)

Wang et al. have simplified the tissue as one-dimensional question [15]; therefore (6a) and (6b) is changed to(30)μ∂2UR∂R2+1R∂UR∂R-URR2=-ρω2u02URμ∂2UZ∂R2+1R∂UZ∂R=-ρω2u02UZ.Based on ∇∙v→=0, the solutions are(31)UR=0UZ=J0ρω2u02μR.Both (29) and (31) do not fit very much in this part; therefore, we suppose (31) as the upper boundary condition (Z=-L0/u0), (29) works as the lower boundary condition (Z = 0), linear interpolation is adopted to calculate part (2), that is,(32)URR,Z=-ZURR,0UzR,Z=UzR,0+ZUzR,0-UzR,-L0u0-L0u0<Z<0.


### 2.6. Parameters

The analytic solution (29) shows that λ affects the amplitude attenuation along z. Langevin et al.’s research on acupuncture demonstrates the displacement is small when z is over 2 cm [10]. Therefore, we suppose the displacement amplitude at z = 2 cm is a small amount compared to the amplitude at z = 0 cm (the plane of the needle tip), that is e-λ(z/u0)≈10%; we can determine λ ≈ 0.25. The standard u0 is 0.5 cm and r0 is 0.02 cm.The biological viscoelasticity is determined by [16] (33)μiω=1+c/2ln⁡1+ω2τ22-ln⁡1+ω2τ12+ictan-1ωτ2-tan-1ωτ11+cln⁡τ2/τ1.where τ1=0.0984,τ2=8454.76,c=0.0351.J 0 ( ρ ω 2 u 0 2 / μ + λ R ) and J1(ρω2u02/μ+λR) will diverge at R→∞ if ρω2u02/μ+λ have an imaginary part. Obviously, λ>0; notice that the imaginary part of ρω2u02/μ is small compared to the real part; for example, ρω2u02/μ = 1.2 − 0.03i when ω=2π and ρω2u02/μ= 76 − 0.38i when ω=16π; therefore we replace ρω2u02/μ+λ with its module ρω2u02/μ+λ to void divergence.


## 2.1. Model

Soft tissues are usually assumed as viscoelastic materials; some are more “elastic” and some are more “viscous”; they exhibit obvious hysteresis, stress relaxation, and other viscoelastic phenomena [14]. In our model, the skin tissue is treated as a semi-infinite viscoelastic body and the needle is treated as a rigid body. As shown in Figure 2, the tissue is divided into two parts. Part (1) represents the tissue from the plane of the tip location to infinity (z ≥0) where tissue displacement is caused by the movement of the needle tip (point A). Part (2) is the needle manipulation region (– L0≤ z < 0), where L0 is the needle insertion depth and the standard L0 is 2 cm in lower limbs.Figure 2
Diagrammatic sketch of lifting-thrusting manipulation.


## 2.2. Equations

The soft tissue is supposed as a homogeneous isotropic incompressible viscoelastic body and fit for Feng's quasilinear viscoelastic theory of small amplitude displacement [15]. For acupuncture manipulation, the needle movement causes small amplitude oscillatory displacements in the low stress range (physiological state) of the tissue. The boundary conditions are specified as simple harmonic functions of time; therefore, the variables can be expressed as(1)uixi,t=u~ieiωtσijxi,t=σ~ijeiωtεijxi,t=ε~ijeiωti,j=1,2,3.where ui is the displacement component which is the distance from the original (equilibrium) position, xi is the coordinate component, t is the time, ω is the angular frequency, σij is the component of stress tensor, and εij is the component of strain tensor.Based on the strain-displacement relation and excluding the effect of volume force, we get the governing equation of motion in linear viscoelastic medium(2)λ∗iω+μ∗iωduj,ji+μ∗iωdui,jj=ρ∂2ui∂t2.where μ∗(iω) and λ∗(iω) are the viscoelastic material parameters andρ is the density. Postulate the tissue is incompressible, that is, uj,ji=0, and putting (1) into (2), we get(3)μ∗iωdu~i,jj=-ρω2u~iRewrite it in tensor form(4)μ∗iω∇2u~→=-ρω2u~→Suppose the tissue around the needle is a semi-infinite space. This is an axisymmetric problem (uθ=0,∂/∂θ=0). Replace μ∗ with μ and u~ with u; then we have (5) in cylindrical coordinate(5)μ∂2ur∂r2+1r∂ur∂r-urr2+∂2ur∂z2=-ρω2urμ∂2uz∂r2+1r∂uz∂r+∂2uz∂z2=-ρω2uz.where ur and uz are the displacement amplitude inr andz direction, respectively. Defining the dimensionless parameters: UZ=uz/u0,UR=ur/u0,Z=z/u0,R=r/u0, where u0 is the amplitude of needle movement. Then (5) is changed to(6a)μ∂2UR∂R2+1R∂UR∂R-URR2+∂2UR∂Z2=-ρω2u02UR(6b)μ∂2UZ∂R2+1R∂UZ∂R+∂2UZ∂Z2=-ρω2u02UZ.Based on the strain-displacement relation(7)εrz=12∂uz∂r+∂ur∂zeiωt=12∂UZ∂R+∂UR∂Zeiωt=ERZeiωt.where ERZ is calculated by explicit difference scheme(8)ERZ=12∂UZ∂R+∂UR∂Z≈12ΔRUZj+1,k-UZj,k+12ΔZURj,k+1-URj,k.And the corresponding stress is(9)τrz=Re2μiωεrz=Reμ1+iμ2·2ERZcos⁡ωt+isinωt=2ERZμ1cos⁡ωt-μ2sin⁡ωt=2ERZμcos⁡ωt+δ.where μ=μ12+μ22, tan⁡δ=μ2/μ1, and δ is the phase difference between stress and strain. The shear stress amplitude TRZ=2μERZ.In a cycle, the energy dissipation of one-unit volume caused by lifting-thrusting is(10)ΔE=∫02π/ωτrzε˙rzdt=4μERZ2∫02π/ω-ωμ1cos⁡ωt-μ2sin⁡ωtsin⁡ωtdt=4πERZ2μ2.The energy dissipation of the whole domain is(11)Ef=∫r0/u0+∞∫-L0/u0+∞ΔEdZdR=∫r0/u0+∞∫-L0/u0+∞4πμ2ERZ2dZdR.Define Emax as the max energy stored(12)Emax=4πERZ2∫0π/ωμ1ωcos⁡ωtsin⁡ωtdt=2πμ1ERZ2.Therefore, the mechanical quality factorQ which reflect the ratio of energy storage to energy dissipation is(13)Q=2πEmaxΔE=μ1μ2.


## 2.3. Boundary Conditions

The motion of the boundary is(14)ur=0,uz=u0cos⁡ωtr=r0,-L0≤z≤0ur=uz=0z→∞uz=0r→∞,wherer0 is the radius of the needle. Then(15)UR=0R→0∪R→∞∪Z→∞UZ=0R→∞∪Z→∞UZ=1R→0∩-L0u0≤Z≤0.


## 2.4. Analytic Solutions of Part (1)

LetUR=F(R)G(Z); (6a) is expressed as(16)FRR+1/RFR-1/R2FF+GZZG+ρω2u02μ=0.Let GZZ/G=λ; then we get(17a)GZZ-λG=0(17b)R2FRR+RFR+ρω2u02μ+λR2-1F=0.Equation (17b) is a 1-order Bessel equation, and the solutions are(18)GZ=A1eλZ+A2e-λZFR=B1J1ρω2u02μ+λR+B2Y1ρω2u02μ+λR,whereJ1,Y1 are 1-order Bessel function and Neumann function, respectively.InputtingUR(Z→∞)=0, we get A1=0. Because J1(0)=0 and Y1(0)=-∞, based on the boundary condition UR(R→0,Z=0)=0, we get B2=0; therefore(19)UR=AJ1ρω2u02μ+λRe-λZ.Because J1(+∞)=0, then UR(R→∞)=0.PutUZ=F(R)G(Z) into (6b), (20)FRR+1/RFRF+ρω2u02μ+GZZG=0.Let GZZ/G=λ; then we get(21a)GZZ-λG=0(21b)FRR+1RFR+ρω2u02μ+λF=0.Equation (21b) is a 0-order Bessel equation; therefore the solutions are(22)GZ=A1eλZ+A2e-λZFR=B1J0ρω2u02μ+λR+B2Y0ρω2u02μ+λR.InputtingUz(Z→∞)=0, we get A1=0. Because J0(0)=1 and Y1(0)=-∞, based on the boundary condition Uz(R→0,Z=0)=1, we get B2=0; therefore(23)UZ=J0ρω2u02μ+λRe-λZ.Suppose the tissue is incompressible, that is,∇∙v→=0; therefore(24)∂UR∂R+URR+∂UZ∂Z=0.Putting (19) and (23) into (24), we get(25)A∂J1ρω2u02/μ+λR∂R+AJ1ρω2u02/μ+λRR-λJ0ρω2u02μ+λR=0.Based on the recurrence relation of Bessel’s function(26)Jv-1x+Jv+1x=2vxJvxJv-1x-Jv+1x=2Jv′x,we get (27)∂J1kR∂R=kJ0kR-J2kR2J1kRR=kJ0kR+J2kR2k=ρω2u02μ+λ.Putting (27) into (25), we get (28)AkJ0kR-J2kR2+J0kR+J2kR2-λJ0kR=0,and then A=λ/k. Therefore(29)UZ=J0ρω2u02μ+λRe-λZUR=λρω2u02/μ+λJ1ρω2u02μ+λRe-λZ.


## 2.5. Calculation Solutions of Part (2)

Wang et al. have simplified the tissue as one-dimensional question [15]; therefore (6a) and (6b) is changed to(30)μ∂2UR∂R2+1R∂UR∂R-URR2=-ρω2u02URμ∂2UZ∂R2+1R∂UZ∂R=-ρω2u02UZ.Based on ∇∙v→=0, the solutions are(31)UR=0UZ=J0ρω2u02μR.Both (29) and (31) do not fit very much in this part; therefore, we suppose (31) as the upper boundary condition (Z=-L0/u0), (29) works as the lower boundary condition (Z = 0), linear interpolation is adopted to calculate part (2), that is,(32)URR,Z=-ZURR,0UzR,Z=UzR,0+ZUzR,0-UzR,-L0u0-L0u0<Z<0.


## 2.6. Parameters

The analytic solution (29) shows that λ affects the amplitude attenuation along z. Langevin et al.’s research on acupuncture demonstrates the displacement is small when z is over 2 cm [10]. Therefore, we suppose the displacement amplitude at z = 2 cm is a small amount compared to the amplitude at z = 0 cm (the plane of the needle tip), that is e-λ(z/u0)≈10%; we can determine λ ≈ 0.25. The standard u0 is 0.5 cm and r0 is 0.02 cm.The biological viscoelasticity is determined by [16] (33)μiω=1+c/2ln⁡1+ω2τ22-ln⁡1+ω2τ12+ictan-1ωτ2-tan-1ωτ11+cln⁡τ2/τ1.where τ1=0.0984,τ2=8454.76,c=0.0351.J 0 ( ρ ω 2 u 0 2 / μ + λ R ) and J1(ρω2u02/μ+λR) will diverge at R→∞ if ρω2u02/μ+λ have an imaginary part. Obviously, λ>0; notice that the imaginary part of ρω2u02/μ is small compared to the real part; for example, ρω2u02/μ = 1.2 − 0.03i when ω=2π and ρω2u02/μ= 76 − 0.38i when ω=16π; therefore we replace ρω2u02/μ+λ with its module ρω2u02/μ+λ to void divergence.


## 3. Results

### 3.1. Boundary(Z=0,Z=-L0/u0) Solutions

Obviously, the upper (Z=-L0/u0) and lower (Z = 0) boundary solutions determine displacement distributions in the skin tissue. The dimensionless tissue displacements (UZ at the upper and lower boundaries and the UR at the lower boundary) at two stimuli frequencies (ω = 2π and ω = 16π) are plotted in Figure 3. The distributions of UZ are very similar at upper and lower boundaries.Figure 3
Amplitudes of tissue displacements at the upper (Z=-L0/u0) and lower (Z = 0) boundaries. (a) ω = 2π; (b) ω = 16π.
(a) (b)Figure4 shows the corresponding absolute values of displacements (UR and UZ) at lower boundary (Z = 0). Additionally, the combined displacement U=UR2+UZ2 is plotted. It shows thatU is nearly the same as UZ and UR is relatively small, especially at high stimulus frequency (ω=16π, in Figure 4(b)). Therefore, we will only studyU afterwards. The influence of stimulus frequency on the displacement (U) is shown in Figure 5. It is clear thatU decreases with the increase ofR, and at higher frequencies, it decreases more rapidly.Figure 4
Amplitudes of tissue displacements atZ = 0. (a) ω = 2π; (b) ω = 16π.
(a) (b)Figure 5
Amplitudes of tissue displacements under different frequencies.The influence of needle movement amplitude (u0) is shown in Figure 6. Figure 6(a) plots the absolute values of tissue displacement (u) at the lower boundary with different u0s (0.25, 0.5, 0.75, and 1.0cm). The range of u is largely dependent on u0, with a nearly linear relationship at R→0. Figure 6(b) shows the dimensionless displacement U=u/u0, which attenuates a little more slowly alongR direction with an increasing u0. Figure 7 demonstrates U in the cases of different needle radius (r0). It shows thatr0 has a minimal effect onU.Figure 6
Amplitudes of tissue displacements under differentu0s at ω = 2π; (a) u. (b) U.
(a) (b)Figure 7
Amplitudes of tissue displacements under differentr0s. (a) ω = 2π; (b) ω = 16π.
(a) (b)


### 3.2. Displacement Amplitude (U) Distribution in the Skin Tissue

Figure8 shows the displacement amplitudeU distribution at ω = 2π and ω = 16π. Comparing Figures 8(a) and 8(b) shows that the tissue displacement around the needle (-L0/u0≤Z≤0) decays along the radial directionR under both low and high frequency stimulus, butU decays less rapidly under low frequency (U ≤ 50% wherer > 0.6 cm) than that under high frequency (U ≤ 50% wherer > 0.07 cm). Moreover,U also decays quickly under the needle tip (Z > 0).Figure 8
Amplitudes of tissue displacements; the black lines represent the displacement atZ=-L0/u0=-4, Z = 0 and Z=L0/u0=4. (a) ω = 2π; (b) ω = 16π.
(a) (b)


### 3.3. Shear Strains (εrz) and Stresses (τrz) under Different Frequencies

Figure9 shows the contour lines of shear strain amplitude in absolute values (ERZ) in the tissue. The needle is located atR= 0, -4≤Z≤0 (shown as white lines in Figure 9). The maximum shear strain occurs at the needle tip, and again, strains decay along the radial direction, with more rapid attenuation under high frequency stimulus. Thus, the low frequency stimulus affects larger area than high frequency stimulus. Take the area where ERZ> 0.2 as the effective influence region Seff, and defineZ andR at ERZ≈0.2 as Zthr and Rthr, respectively; then Seff =Zthr×Rthr≈1×0.94 atω = 2π and Seff ≈ 1×0.08 atω = 16π, shown as the black rectangles in Figure 9. It turns out that the stimulus frequency has little effect on Zthr but has significant influence on Rthr. Figure 10 shows the contour lines of dimensionless shear stress amplitude in absolute values (TRZ); it shows the same tendency as ERZ.Figure 9
Contour lines ofERZ in the tissue; the white lines at r→0∩-4≤Z≤0 represent the needle; the black rectangles represent the approximate area where ERZ > 0.2. (a)ω = 2π; (b)ω = 16π.
(a) (b)Figure 10
Contour lines ofTRZ in the tissue, the white lines at R→0∩-4≤Z≤0 represent the needle. (a)ω = 2π; (b)ω = 16π.
(a) (b)Table1 lists Rthrs and Seffs under different frequencies. Figure 11 shows that there is an approximately linear relationship between Rthrs and the reciprocal of the stimulus frequency 1/ωs. The same relationship also exists between Seffs and 1/ωs.Table 1
R t h rs and Seffs under different angular frequencies.
ω 2π 4π 6π 8π 10π 12π 14π 16π R t h r 0.94 0.56 0.36 0.26 0.20 0.16 0.12 0.08  S eff 0.5×0.94 0.5×0.56 0.5×0.36 0.5×0.26 0.5×0.20 0.5×0.16 0.5×0.12 0.5×0.08Figure 11
Relationship betweenRthrs and Seffs on one side and 1/ωs on the other. Diamonds represent calculating results and blue lines represent the fitting curve.
(a) (b)


### 3.4. Energy Dissipation (ΔE) and Mechanical Quality Factor (Q) under Different Frequencies

Figure12 shows ΔE andQ under different frequencies. Figure 12(a) shows that ΔE is larger and more intensified around the needle at high frequencies than that at low frequencies. For the mechanical quality factorQ, a nearly linear relationship between Qs and the stimulus frequencies can be found (shown in Figure 12(b)). This means that energy storage ratio increases with the increase of stimulus frequency, and the tissue demonstrates more “elastic” at high frequencies.Figure 12
ΔE andQ under different frequencies. (a) ΔE at Z = 0; (b) Q.
(a) (b)


## 3.1. Boundary(Z=0,Z=-L0/u0) Solutions

Obviously, the upper (Z=-L0/u0) and lower (Z = 0) boundary solutions determine displacement distributions in the skin tissue. The dimensionless tissue displacements (UZ at the upper and lower boundaries and the UR at the lower boundary) at two stimuli frequencies (ω = 2π and ω = 16π) are plotted in Figure 3. The distributions of UZ are very similar at upper and lower boundaries.Figure 3
Amplitudes of tissue displacements at the upper (Z=-L0/u0) and lower (Z = 0) boundaries. (a) ω = 2π; (b) ω = 16π.
(a) (b)Figure4 shows the corresponding absolute values of displacements (UR and UZ) at lower boundary (Z = 0). Additionally, the combined displacement U=UR2+UZ2 is plotted. It shows thatU is nearly the same as UZ and UR is relatively small, especially at high stimulus frequency (ω=16π, in Figure 4(b)). Therefore, we will only studyU afterwards. The influence of stimulus frequency on the displacement (U) is shown in Figure 5. It is clear thatU decreases with the increase ofR, and at higher frequencies, it decreases more rapidly.Figure 4
Amplitudes of tissue displacements atZ = 0. (a) ω = 2π; (b) ω = 16π.
(a) (b)Figure 5
Amplitudes of tissue displacements under different frequencies.The influence of needle movement amplitude (u0) is shown in Figure 6. Figure 6(a) plots the absolute values of tissue displacement (u) at the lower boundary with different u0s (0.25, 0.5, 0.75, and 1.0cm). The range of u is largely dependent on u0, with a nearly linear relationship at R→0. Figure 6(b) shows the dimensionless displacement U=u/u0, which attenuates a little more slowly alongR direction with an increasing u0. Figure 7 demonstrates U in the cases of different needle radius (r0). It shows thatr0 has a minimal effect onU.Figure 6
Amplitudes of tissue displacements under differentu0s at ω = 2π; (a) u. (b) U.
(a) (b)Figure 7
Amplitudes of tissue displacements under differentr0s. (a) ω = 2π; (b) ω = 16π.
(a) (b)


## 3.2. Displacement Amplitude (U) Distribution in the Skin Tissue

Figure8 shows the displacement amplitudeU distribution at ω = 2π and ω = 16π. Comparing Figures 8(a) and 8(b) shows that the tissue displacement around the needle (-L0/u0≤Z≤0) decays along the radial directionR under both low and high frequency stimulus, butU decays less rapidly under low frequency (U ≤ 50% wherer > 0.6 cm) than that under high frequency (U ≤ 50% wherer > 0.07 cm). Moreover,U also decays quickly under the needle tip (Z > 0).Figure 8
Amplitudes of tissue displacements; the black lines represent the displacement atZ=-L0/u0=-4, Z = 0 and Z=L0/u0=4. (a) ω = 2π; (b) ω = 16π.
(a) (b)


## 3.3. Shear Strains (εrz) and Stresses (τrz) under Different Frequencies

Figure9 shows the contour lines of shear strain amplitude in absolute values (ERZ) in the tissue. The needle is located atR= 0, -4≤Z≤0 (shown as white lines in Figure 9). The maximum shear strain occurs at the needle tip, and again, strains decay along the radial direction, with more rapid attenuation under high frequency stimulus. Thus, the low frequency stimulus affects larger area than high frequency stimulus. Take the area where ERZ> 0.2 as the effective influence region Seff, and defineZ andR at ERZ≈0.2 as Zthr and Rthr, respectively; then Seff =Zthr×Rthr≈1×0.94 atω = 2π and Seff ≈ 1×0.08 atω = 16π, shown as the black rectangles in Figure 9. It turns out that the stimulus frequency has little effect on Zthr but has significant influence on Rthr. Figure 10 shows the contour lines of dimensionless shear stress amplitude in absolute values (TRZ); it shows the same tendency as ERZ.Figure 9
Contour lines ofERZ in the tissue; the white lines at r→0∩-4≤Z≤0 represent the needle; the black rectangles represent the approximate area where ERZ > 0.2. (a)ω = 2π; (b)ω = 16π.
(a) (b)Figure 10
Contour lines ofTRZ in the tissue, the white lines at R→0∩-4≤Z≤0 represent the needle. (a)ω = 2π; (b)ω = 16π.
(a) (b)Table1 lists Rthrs and Seffs under different frequencies. Figure 11 shows that there is an approximately linear relationship between Rthrs and the reciprocal of the stimulus frequency 1/ωs. The same relationship also exists between Seffs and 1/ωs.Table 1
R t h rs and Seffs under different angular frequencies.
ω 2π 4π 6π 8π 10π 12π 14π 16π R t h r 0.94 0.56 0.36 0.26 0.20 0.16 0.12 0.08  S eff 0.5×0.94 0.5×0.56 0.5×0.36 0.5×0.26 0.5×0.20 0.5×0.16 0.5×0.12 0.5×0.08Figure 11
Relationship betweenRthrs and Seffs on one side and 1/ωs on the other. Diamonds represent calculating results and blue lines represent the fitting curve.
(a) (b)


## 3.4. Energy Dissipation (ΔE) and Mechanical Quality Factor (Q) under Different Frequencies

Figure12 shows ΔE andQ under different frequencies. Figure 12(a) shows that ΔE is larger and more intensified around the needle at high frequencies than that at low frequencies. For the mechanical quality factorQ, a nearly linear relationship between Qs and the stimulus frequencies can be found (shown in Figure 12(b)). This means that energy storage ratio increases with the increase of stimulus frequency, and the tissue demonstrates more “elastic” at high frequencies.Figure 12
ΔE andQ under different frequencies. (a) ΔE at Z = 0; (b) Q.
(a) (b)


## 4. Discussion

Langevin et al. claimed that acupuncture manipulation would lead to tissue displacement and generate a mechanical stress field [10]. This phenomenon is well demonstrated in our simulation. Moreover, the local displacement, strain, stress, and energy near the acupoint (or the needle) have been quantitively analyzed. The characteristics of the tissue displacement are concluded as (1) it is larger at low frequencies (see Figure 8) given other conditions the same (the needle movementu0 and the material properties, etc.), and (2), the displacement will decay more slowly and extend relatively further from the needle at low frequencies (see Figure 5).Main characteristics in the mechanical process of acupuncture are the stimulus amplitude (u0) and the angular frequency (ω). However, there is a lack of consistency of these characteristics. Wang XM pointed out that mechanical information can only influence nearer region of the acupoint [17]. According to our simulation, the resultant tissue displacement is nearly proportional tou0 (see Figure 6(a)). It is also shown that the displacement, strain, and stress under low frequency stimulus decay more slowly than those under high frequency (see Figures 8–10). These results indicate that the mechanical information (displacement, strain, and stress) can propagate further under low frequency stimulus. Figure 13 shows the characteristic curves of different acupuncturists, where significant differences in u0 and ω can be found, but u0 is larger at high frequencies; thus u0/ω tends to be constant [18]. Figure 11(a) shows that there is nearly an inverse relationship between Rthrs andωs, or Rthr=k/ω, where k is a constant. Given R=r/u0, then the radius of the effective influence area rthr=u0×Rthr=u0×k/ω=k(u0/ω). Thus a constant u0/ω is equivalent to a constant rthr; therefore, to get good acupuncture effect, a sufficient area of tissue should be deformed, which may be an alternative interpolation of acupuncture sensations in the mechanical viewpoint. However, our model is not able to figure out what “sufficient area” should be, which requires further experimental research. But suggestions on achieving mechanical sensations can be made. Since needle movement is a damage process to the tissue and slip between the needle movement and tissue movement may occur at large u0 (we suppose there is no slip because of small u0), thus u0 is limited. Therefore, to achieve better effect (larger effective influence region), low frequency stimulus is a good choice. It was reported that the frequency of acupuncture sensations is concentrated at low frequency [18]. Ding et al. reported the mean frequency to achieve acupuncture sensations by experienced acupuncturists is 1.2 Hz [19].Figure 13
Characteristic curves of five professional acupuncturists.Figure7 shows that needle diameter has no effect on displacement. Li et al. reported needling of different diameters induced the consistent change of fascial connective tissues and collagenous fiber arrangement [20]. Figures 9 and 10 show that strain and stress decay quickly under the needle tip. That means the needle must be inserted into the right depth (about 2 cm in lower limbs) and position (π/ω cm around the acupoint in our simulation) to induce enough strain and stress. A survey showed that acupuncturists think deep insertion was easier to induce acupuncture sensations than shallow insertion [21].Because of the limitation of experimental means, it is impossible to measure tissue strain and stress in vivo, tissue displacement, or the stress on the needle are exploited to infer the tissue strain [22]. But there is a complicated relationship between strain and displacement; it is controversial to infer stress field without considering needle frequency. Our model predicted the relationship between the main deformation areas and frequencies (see Figure 11). But there is no agreement on the threshold of ERZ to achieve acupuncture effect and tissue is more complicated (anisotropy, nonlinearity) than our model [23]; further experiments are required to testify the results.Figure12(a) shows that energy dissipation is relatively small under low frequency compared with that under high frequency; Figure 10(b) shows that there is a linear relationship between mechanical quality factors and frequencies. Obviously, a high energy is required to maintain a high frequency manipulation.


## 5. Conclusion

In this study, we have illustrated that lifting-thrusting needle will lead to tissue displacements and generate a mechanical stress field. The displacement correlated with needle frequency and decayed more slowly under low frequency than that under high frequency. Strain and stress correlated with frequency too, and there exists an inverse relationship between the main deformation areas and frequencies. Energy dissipation is also affected by frequency. This is a preliminary study of mechanical effects of acupuncture manipulation, isotropy, and quasi-linear viscoelastic theory adopted to simplify soft tissue, while Fox et al. observed anisotropic tissue motion induced by acupuncture [24]. Our results may reflect local effects around the needle; anisotropy must be considered if we study the mechanical information propagation in the tissue.



---

*Source: 5475426-2019-04-08.xml*