### Robot

A robotic 6-degree-of-freedom setup (KUKA KR 60–3 robot, Augsburg, Germany; reproducibility: ±0.06 mm) including a universal force/torque sensor (ATI UFS: Theta SI1000-120; resolution: 0.25 N and 0.025 Nm) was used to perform axially loaded knee flexion.

### Specimen preparation

Four (two left and two right) fresh frozen, healthy, skinless sheep knee joints were obtained post mortem and directly stored at -20°C. Approximate body weights of the four sheep were 30 kg, 30 kg, 30 kg and 35 kg. Prior to testing, the joints were thawed for about 16 hours at room temperature wrapped in a cloth soaked with saline solution. The femur and tibia were resected 20 cm proximal and distal to the joint gap and embedded in a two-component resin (RenCast^{©} FC 53 isocyanate/FC 53 polyol, Gößl & Pfaff GmbH, Karlskron, Germany). The embedded bones were fixed within aluminum cylinders using radial screws around the circumference. For the testing procedure, the tibial cylinder was fixed to a column stand which was connected to the floor while the femoral cylinder was attached to the robotic system. During experiments, the specimens were wrapped in thin polyethylene film to preserve the specimens from drying (Figure
1). After the experiments, the knee joints were dissected to visually control the joint structures.

### Axes definitions

The femoral epicondyles of the knee joint were palpated and marked with a black felt pen. The initial origin of the robotic coordinate systems of the base (tibia, fixed) and tool (femur, moving) were placed in the midpoint between the lateral and medial epicondyle. For the tibial and femoral local coordinate systems, the x-axes were defined from medial to lateral, the y-axes from posterior to anterior and the longitudinal z-axes from distal to proximal. The x-axis for each of the two coordinate systems was a flexion axis specified perpendicular to the respective femoral and tibial longitudinal z-axis, approximately coinciding with the epicondylar axis (Figure
1). Tibiofemoral movement was recorded in the tibia (base) coordinate system.

### Recording of passive knee flexion

Similar to the previous study
[3], passive knee flexion was recorded with seven different increments (0.1, 0.2, 0.3, 0.4, 0.5, 0.75 and 1°) for one specimen in order to determine the optimal angular resolution for the knee flexion motion. Passive knee flexion of another three specimens was recorded with the optimal increment. The optimal increment was determined as a minimal increment with which the desired flexion velocity could be performed by the robot.

Passive knee flexion was characterized in a way that it follows the path of minimal resistance. It has been shown that this path of passive flexion was unique for each knee joint and therefore describes the individual unloaded motion of the joint
[12, 13]. Passive knee flexion was recorded with the robotic system using a combination of angle controlled and force/torque controlled motion
[14]. While the flexion movement around the local femur flexion medio-lateral axis was carried out in angle controlled mode, the remaining five degree of freedom were force/torque controlled at the same time. All forces and torques were specified to be zero (acceptable tolerance for the residual magnitude less than 2 N and 0.2 Nm respectively), except for the axial force (along the tibia longitudinal axis), where 10 N were applied to guarantee a contact of the two joint surfaces
[15]. The knee flexion motion was started by approximately 50° of the knee flexion angle which corresponds to the minimal knee flexion angle during the heel-strike recorded by Taylor et al.
[11], then the knee was flexed for 30° with the specified increment (until the flexion angle of 80°). Each point of the passive flexion path was recorded after all force/torque conditions were met. Thus, for instance, after recording passive knee flexion with an increment of 0.5°, the path of the passive knee flexion consisted of 61 points (an initial point plus 60 points after rotations with an increment of 0.5°).

In order to optimize the robot movement around the knee flexion axis, the origins of the tibial (base) and femoral (tool) coordinate systems were newly defined in the knee rotation center. Therefore, passive knee flexion was performed twice: the first time using a manually determined knee rotation center and the second time using the calculated knee rotation center. Thereafter, Matlab software (Matlab R2013a, The Math-Works, Inc., Natick, MA, USA) was used for the data processing described below.

The knee rotation center was calculated as a minimal amplitude point
[16] during tibiofemoral movement of the first recorded passive knee flexion. Ehrig et al. showed that the minimal amplitude method was less prone to noise
[17] than the least-squares algorithm which was used for the same purpose in the previous study
[3]. Thus, the robotic flexion axis was corrected to the calculated knee rotation center which corresponded to that optimized axis and the passive knee flexion was recorded again.

### Preparation of dynamics data for a robot-assisted knee flexion simulation

The curves of averaged sheep’s knee flexion and contact axial force were taken from the work of Taylor et al.
[11] to simulate an in vivo load on an ovine knee specimen during the flexion. The data of both curves from heel-strike to heel-strike (100% of a complete gait cycle) were corrected at the start and end of the gait cycle in order to avoid a pulse movement of the robot during a juncture of two adjacent gait cycles. To perform this correction, the data of three gait cycles were connected in one curve and processed with a low pass filter and then the corrected data of the second gait cycle were used to interpolate flexion and axial force curves for the robot.The flexion angle data were interpolated for all of the seven increments. An example of the incremental interpolation of the flexion angle data was shown in Figure
2A. If two adjacent, interpolated points at local minima or maxima had the same flexion angle, the midpoint was taken instead of both points (Figure
2A, black arrow) in order to avoid zero angular velocity between the two points.The absolute angular velocity to rotate the femur around the x-axis by means of the robot (Figure
1) was calculated as a differential of the knee flexion angle data for each of the seven increments. To undergo the differentiation, the end value of the knee flexion angle data was added to the flexion angle data set as its first element in order to avoid shortening of the data set length after the differentiation. Then, the data set of the absolute angular velocity values was normalized to a maximal value from this data set. Thus, to calculate an angular velocity profile for a gait cycle with a desired maximal angular velocity, the set of the normalized angular velocity data had to be multiplied with the desired maximal value.

The axial contact force values in body weight (BW) for each of the seven increments were interpolated from the knee contact force values published by Taylor et al.
[11] to corresponding gait cycle values of the flexion angle curves (see the example shown in Figure
2B). Thus, to calculate a set of axial forces for a gait cycle with a desired body weight in Newton (N), the set of the axial contact force data in BW had to be multiplied with the sheep’s body weight in N.

### In vivo loaded knee flexion simulation

The flexion angle, angular velocity and axial force data sets were vectors of the same length. These three vectors of data comprised dynamics information for the simulation of knee flexion over one gait cycle. The loaded knee flexion path of the robot movement consisted of points (spatial positions of the robot) from the previously recorded passive flexion path. In order to reproduce the knee flexion during the gait cycle, we substituted the values of the interpolated flexion angle data for indexes of the corresponding points of the passive flexion path. Thus, the loaded knee flexion path was constructed from points of the passive flexion path to simulate the sheep’s knee flexion.

The robot performed the loaded knee flexion path in a superimposed force torque control mode. Only the axial force was actively simulated (controlled) by the robot. Medio-lateral and posterior-anterior forces resulted as a backlash of knee joint structures. The simulated body weights were 300 N, 300 N, 300 N, and 350 N corresponding to the sheep’s body weights. The in vivo loads were simulated for 200 gait cycles on each of the four knee specimens. The 60 initial gait cycles were considered as pre-conditional cycles, while the last 140 gait cycles were used for further data analysis.

### Experiment with flexion velocity over 10°/s limit

As a precaution of our occupational health and safety department, the angular velocity of the robot was limited to 10°/s. This allowed work with the robot without a protective fence. In order to show the reaction of the robot control to target angular velocities above this 10°/s limit, the knee flexion was performed on one specimen with three different flexion velocity profiles. Three flexion velocity profiles were calculated for maximal flexion velocity of 10, 15 and 20°/s.

### Comparison of simulated forces to in vivo forces

The dynamics data of simulated knee flexion were recorded with a sampling time of 24 ms. In order to calculate the residual error between the simulated force data and the in vivo force data published by Taylor et al.
[11], the in vivo force data were interpolated for each percent of the gait cycle while the simulated force data were averaged for a category of 1% during the gait cycle. Thus, the residual error *ϵ* in BW units was calculated as follows:

\u03f5=\sqrt{\frac{1}{m\times n}{\displaystyle \sum _{j=1}^{m=4}{\displaystyle \sum _{i=1}^{n=140}{\left({F}_{ij}-{F}_{\mathit{invivo}}\right)}^{2}}}}\text{,}

where m denoted four specimens, n denoted 140 gait cycles. *F*
_{
ij
} was the simulated force vector and *F*
_{
invivo
} was the in vivo force vector corresponding to Taylor et al.
[11]. Both force vectors consisted of 100 elements in body weight units.

The total residual error *ϵ*
_{
total
} in BW units was calculated as follows:

{\u03f5}_{\mathit{total}}=\sqrt{\frac{1}{p\times m\times n}{\displaystyle \sum _{k=1}^{p=100}{\displaystyle \sum _{j=1}^{m=4}{\displaystyle \sum _{i=1}^{n=140}{\left({F}_{kij}-{F}_{k\_\mathit{invivo}}\right)}^{2}}}}}\text{,}

where *p* denoted 100 elements of a force vector. *F*
_{
kij
} was an element of the simulated force and *F*
_{
k _ invivo
} was an element of the in vivo force.

### Standard deviation of the simulated forces

The simulated force data were averaged for a category of 1% during the gait cycle. Thus, the gait cycle consisted of 100 elements. For each percent of the gait cycle, a standard deviation of the forces was calculated in body weight units. The total standard deviation *σ*
_{
total
} in BW units was calculated as follows:

{\sigma}_{\mathit{total}}=\sqrt{\frac{1}{p\times m\times n-1}{\displaystyle \sum _{k=1}^{p=100}{\displaystyle \sum _{i=1}^{m\times n=560}{\left({F}_{ki}-{\overline{F}}_{k}\right)}^{2}}\text{,}}}

where
{\overline{F}}_{k}=\frac{{\displaystyle \sum _{i=1}^{m\times n=560}{F}_{ki}}}{m\times n}, *m* denoted four specimens, *n* denoted 140 gait cycles and *p* denoted 100 elements of a force vector.