As the primary body weight support mechanism, the musculoskeletal system is sensitive to the impact loading of the joints and the bones during walking and galloping, the analysis of the mechanical properties of the bones and the soft tissues can provide a reference to the design of bionic mechanisms. Based on the anatomical data, we construct a mathematical model for the musculoskeletal system of a domestic cat.

## 2.1 Musculoskeletal system modeling:

In the anatomical data, the truck of the cat is 280 mm long, the shoulder is 60/70 mm wide, and the body weight is 1.4 kg. The other measured data is given in Table 1.

Table 1

Parameters of the limb bones and muscles

skeleton and muscle | length/mm | diameter /mm |

scapula | 44 | -- |

humerus | 69 | 6.8 |

radius | 58 | 4 |

metacarpal | 22 | 3.1 |

femur | 80 | 6.9 |

tibia | 82 | 6 |

metatarsal | 35 | 3 |

Posterior tibia muscle group | 90 | -- |

Lateral femoral muscle group | 85 | -- |

Posterior femoral group muscle group | 90 | -- |

The distribution of the tendon attachment data is also obtained from the anatomical measurements, where the muscles are distributed along various joints and bones, and the muscle forces are correlated to the motions of different bones and joints to drive the bones and joints swinging, transforming the muscle contractions into the rotation around the joint of the bones. Under the parameters shown in Table 1, the location of the attachment point of muscle tendon along the bone can be utilized to construct the kinematical model of the musculoskeletal system, while the bones and joints are simplified according to their functions of the body activity in this paper. Furthermore, the coordinate transformation method is employed for the kinematical modeling of the musculoskeletal system, as shown in the following equations:

$${T_{hind\lim b}}={T_{hip}}{T_{knee}}{T_{ankle}}{T_{toe}}$$

1

Eq. (1) is the kinematical equations of hind limbs, containing the relationship of the coordinate transformation from the toe end to the hip joint, where **T****hip**, **T****knee**, **T****ankle**, **T****toe** are the kinematical transform equations of the latter joint relative to the former joint. The kinematic equations of the other limbs are similar.

**Mradius**: Torque acting on radius; **Mhumerus**: Torque acting on humerus; **Mmetacarpal**: Torque acting on metacarpal; **Mscapula**: Torque acting on scapula; **Mtibia**:Torque acting on tibia; **Mfemur**:Torque acting on femur: **Mtarsal**:Torque acting on tarsal

Using the kinematical model of the musculoskeletal system (shown in Fig. 1b) analyzed above and by measuring the running pattern[35, 36] of the cat shown in Fig. 1a, the change in each joint angle can be obtained. The change of the muscle length and contraction velocity can be solved analytically as well. The muscle forces can be calculated by solving Eq. (2):

$${F_m}={F_a}+{F_p}={F_0}({f_1}{f_2}a(t)+{f_3})$$

2

In Eq. (2), the distance change between the muscle tendon attachment points in Fig. 1b varying with time is known, **F****m** represents the size of the muscle force represented by the single connection lines shown in Fig. 1b could be obtained., **F****a** and **F****p** represent the forces of the flexor muscle and extensor muscle, respectively. **f****1**, **f****2**, **f****3** represent the exponential functions where are the inverse tangent function of the change of the muscle length and contraction velocity, respectively.

Finally, the muscle force acting at each joint within the running time of 0.3 s is calculated. The muscle forces of the hind limbs are generally larger than those of the fore limbs, which shows that the hind limbs are the major contributors to body support during running.

## 2.2 Determination of biomaterial properties

In Fig. 1c, the torque acting on different joints is obtained by multiplying the muscle force (calculated in Eq. (2)) and the moment arm, where the moment arm is the minimum distance between the connection line (shown in Fig. 1b) and the rotation center of the joint[32]. As described above, hind limbs exert the major power for high speed movement[33], and tibias bridge the preceding and the following for the transmission of movement. The mechanism of tibias can bear a certain bending moment and torque moment, can also withstand the radial impact pressure. Therefore, in this paper we select the tibia to study the dynamic response of bones of the domestic cat under impact loading.

The tibia used in this study is 8.3 cm long, and its lateral ankle is 2.2 cm long and weighs 1.41 g, with a density of about 298.7 kg/m3. The truncated tibia is 4.1 cm long and weighs 1 g with a density of about 1940.9 kg/m3. In order to precisely calculate the dynamic response, a nano-indentation technique is employed to measure Young's moduli of the tibia. Figure 2 depicts the displacement-loading curves of the tibia material, where the abscissa represents the pressed depth, and the ordinate represents the press-in force. Figure 2a shows the result for the outside surface of tibia, and Fig. 2b for the inside of the tibia and the medial condyle. There are eight groups of reasonable measured values in Fig. 2a, six in Fig. 2b. We feed the measuring results into the O&P method[34] to obtain multi-group elastic moduli of tibia material, which are averaged to get the elastic moduli of the truncated tibia and the medial condyle as 7.45Gpa and 1.485Gpa respectively.

## 2.3 Numerical model of the dynamic response process

To investigate the transient response of bones in this paper, the contact between the tibia and the tarsal is simplified: the tibia is simplified into an one-dimensional plane rod, and the tarsal into a rigid surface, as shown in Fig. 3a, where under the action of the joint's moment **M**, the femur exerts a pressure on the tibia, and a brief contact-impact occurs between them to generate the dynamic response.

To handle different elastic moduli between the cancellous bone and the compact bone, as shown in Fig. 3b, the tibia is simplified into the collision rod of different materials, where the lower end of the rod represents the condyle. The cancellous bone is relatively abundant. The elastic moduli is *E**1*. In the middle of tibia the elastic moduli is *E**2*, and there is more compact bone.

To solve the nonlinear problem of the dynamic response, the substructure technique for dynamics is employed in this paper. The method can reduce the modal order of the finite element method, and ensure the accuracy of the calculation. Based on the substructure method, we divide the collision rod into *n* substructure rods, each containing *m* two-node link elements, as shown in Fig. 3c, where *C* represents a contact point, **a**, **V** is the acceleration and speed of the rod, respectively, *u* the physical displacement of the node. The motion equation of the substructure **S****(k)** can be expressed as:

$$\left[ {\begin{array}{*{20}{c}} {M_{{ii}}^{{\left( s \right)}}}&{M_{{ib}}^{{\left( s \right)}}} \\ {M_{{bi}}^{{\left( s \right)}}}&{M_{{bb}}^{{\left( s \right)}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\ddot {u}_{i}^{{\left( s \right)}}} \\ {\ddot {u}_{b}^{{\left( s \right)}}} \end{array}} \right\}+\left[ {\begin{array}{*{20}{c}} {K_{{ii}}^{{\left( s \right)}}}&{K_{{ib}}^{{\left( s \right)}}} \\ {K_{{bi}}^{{\left( s \right)}}}&{K_{{bb}}^{{\left( s \right)}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_{i}^{{\left( s \right)}}} \\ {u_{b}^{{\left( s \right)}}} \end{array}} \right\}=\left\{ {\begin{array}{*{20}{c}} {F_{i}^{{\left( s \right)}}} \\ {F_{b}^{{\left( s \right)}}+R_{b}^{{\left( s \right)}}} \end{array}} \right\}$$

3

Namely: (4)

The mass matrix of Eq. (4) is a coordinated mass matrix, where **F** is the external load, and **R** the interfacial force. The physics displacement vector **u****(s)** is divided into two parts: the internal node displacement **u****i****(s)** and the interface node displacement **u****b****(s)**. By calculating the transforming relationship between the physical displacement and the modal displacement, the modal matrix can be obtained, and Eq. (4) can be solved based on the known parameter matrices. Meanwhile, when the rod vertically impacts the ground, we assume no invasion between the rod and the ground, and no geometric dispersion along the rod. As shown in Fig. 3c, during the solution process, the rod length *L* is divided into about 300 units, with every four nodes making a substructure unit. The convergence of the substructure method has been proved.