Ambient Pressure Variation

The ambient pressure variation is obtained from Eq. 1 for the geometrical parameters specified in the modelling section. For the geometry considered, at a load of 13 N the pressure in the synovial fluid goes well below its vapour pressure required to initiate cavitation and is taken as the load on the synovial fluid for the baseline simulation. The evolution of the ambient pressure in the joint for the baseline simulation is shown in Fig. 2a. Clearly, almost all the changes in the ambient pressure occur within the first 0.01 seconds, making the model for the ambient pressure consistent with the previous experimental observations of Unsworth et al.2. Furthermore, the ambient pressure grows increasingly fast initially, before gradually saturating at the zero load pressure when the joint separation becomes large. This has important consequences for the dynamics of the cavitation bubble.

Figure 2 (a) Variation of the ambient pressure in the joint and the excess pressure of the bubble inside the joint from the solution of Eqs 1 and 5, respectively for an acceleration a of 72 m/s2. The inset shows the eccentricity of the joint at 1 ms, the time up to which an excellent match was obtained with experiments. (b) Solution of the Rayleigh-Plesset equation (Eq. 5) for the ambient pressure obtained from Eq. 1 for an acceleration a of 72 m/s2. (i) Evolution of the normalized bubble wall radius for the simulated parameters. (ii) Evolution of the normalized bubble wall velocity. Full size image

Dynamics of the Cavitation Bubble

For resolving the dynamics of the cavitation bubble, the Rayleigh-Plesset equation was coupled to the obtained ambient pressure variation and solved for a bubble with an initial radius of 200 μm. The variation of the bubble radius and bubble wall velocity with time are shown in Fig. 2b. As expected, in response to the ambient pressure variation, the bubble radius decreases, indicating bubble collapse. The negative values of bubble wall velocity also reflect the collapse, although it is evident that the rate of change of this velocity decreases and changes direction. This occurs when the pressure inside the collapsing bubble builds up faster than the rate at which the ambient pressure increases (Fig. 2a), causing the bubble collapse to decelerate. It is worth noting that only the collapsing phase of the bubble (i.e. negative wall velocities) is simulated and that the bubble persists following a partial collapse.

Acoustics

The collapsing bubble produces acoustic pressure waves. The magnitude of these waves is obtained by solving the acoustic pressure equation (Eq. 7) for the collapsing bubble assuming the measurement point to be 0.01 meters from the centre of the bubble (approximate separation of the microphone from the MCP joint in the experiments). The acoustic waveform obtained from the experiments displays a characteristic dominant peak followed by a series of damped oscillations. A peak magnitude of 83 dB and a dominant frequency of 129 Hz are obtained for the acoustic pressure waves from the baseline simulation, which is comparable to the experimental reports in literature3.

Comparison with Experiments

Figure 3 depicts a comparison between the computed normalized acoustic pressure waveform and the waveforms obtained from experiments. The simulated results are seen to correlate well with the experimental measurements, exhibiting an initial dominant pressure peak around 0.1 ms (which as in Fig. 2a corresponds to an instantaneous eccentricity of 3.38 mm) after the initiation of the cracking followed by a series of damped pressure waves. The subsequent peaks after the initial dominant peak are however under predicted in the simulations. Nonetheless, the audio signature of the simulations compares very well with that of the experiments as the audio signature is dominated by the characteristics of the initial peak [See supplementary video]. Hence, given the variability in the acoustic pressure waveforms across subjects and even for the same subject at different instances (Fig. 3 green and black curves), the relatively simple model proposed here appears to be capable of capturing the essential features of the sounds accompanying knuckle cracking.

Figure 3 (a) Normalized acoustic pressure predicted by the model (blue curve) and measured experimentally (other curves). All curves exhibit an initial peak followed by a series of damped pressure oscillations. The inset illustrates the reasonably good agreement in the behaviour of the initial dominant peak which is responsible for the cracking sound. (b) Normalized amplitude-frequency spectrum for the simulation and the experiments. The inset shows that the dominant frequency obtained from the simulations is consistent with the frequencies obtained from the experiments. Full size image

To further characterize these pressure waveforms, we performed a Fast Fourier Transform (FFT) and compared the frequency spectra of the simulations to those from experiments. Figure 3b demonstrates that the frequency at which the largest amplitude is observed is quite similar in the simulations and the experiments. Since the acoustic signature is dominated by the frequency (and its harmonics) having the highest amplitude, the acoustic signature of the simulations, as previously stated, is well correlated to those measured experimentally.

The good correlation between the simulations and experiments also suggests that bubble collapse is a plausible reason for the cracking sound. This result is in agreement with the conclusions of Unsworth et al.2. The competing theory of tribonucleation-mediated bubble inception as the source of the sound stems from the observation that stable bubbles remain in the synovial fluid after an audible sound is heard4. Interestingly, our results show that a partial collapse of the cavitation bubble leads to physiologically consistent acoustics and leaves behind a stable micro-bubble, suggesting that the cavitation theory of knuckle cracking acoustics could indeed explain the experimental observation of stable micro-bubbles. The final micro-bubble in our simulations remains unchanged infinitely long as we assumed an impervious bubble wall. However, future numerical studies may model bubble inception and employ a more realistic bubble wall permeability to confirm and study the inception, terminal and long term dynamics of these bubbles. If the the role of bubble inception on acoustics is tested and the presence of these terminal micro-bubbles after a cavitation bubble collapse is confirmed, all the experimental observations in the literature can be reconciled and thus the disagreements surrounding the source of the sound can be resolved.

Sensitivity Analysis

To analyse the robustness of the model to the real-life variations of parameters across different knuckles, we performed a sensitivity analysis of the acoustic signal to the various parameters in the model including the joint acceleration, initial cavitation bubble radius, fluid viscosity, and joint spacing. To this end, we began by performing an order-of-magnitude analysis on the non-dimensionalized equations which revealed that the pressure term in the Rayleigh-Plesset equation had the strongest influence on the dynamics of the system (see Methods section). Hence, the parameters that influence the pressure term such as the joint acceleration, the load acting on the synovial fluid and the joint geometry are expected to have the largest effect on the acoustic signature.

Figure 4a illustrates the dependence of the acoustic spectra on joint acceleration (a), the initial bubble radius (R 0 ) and the load (w). The acceleration was varied from 36 m/s2 to 144 m/s2, accounting for the estimated experimental uncertainties in the spatial (±0.3 mm) and temporal (±0.002 s) resolution of the joint position5. An increase in the joint acceleration is seen to increase the relative magnitude of the high frequency components in the spectra, without affecting the dominant frequency. The frequency spectra are also relatively invariant to changes of ±100 μm in initial bubble radius. However, the acoustic signature is sensitive to changes in the applied load, with a change in the load from 10 N to 14 N changing the dominant frequency from 64 Hz to 193 Hz. More interestingly, the parametric study on the effect of the load showed that for every geometry considered, there was a limiting value of the load beyond which the model failed to converge. This limiting value is a function of the geometry (both the radius of curvatures and the half angular extent of the joint) and for the cases we considered, we found the maximum limiting value to be around 15 N or about 15% of the experimentally measured external load at which the joint cracks1,2. Acoustic signatures close to the experiments were obtained for a load of ~13 N, which is slightly less than this maximum limiting value. We hypothesize that this threshold provides a ballpark estimate for the hypothesis advanced by Unsworth et al.2 that only a fraction of the applied external load acts on the synovial fluid, with the rest being taken up by the surrounding tissue.

Figure 4 Sensitivity analysis of the model. (a) The model is sensitive to changes in the load but is relatively invariant to changes in acceleration and bubble radius. (b) The model gives physiologically consistent results but is also sensitive to changes in the geometry as shown by varying the geometry through half a standard deviation (σ m = 1.08, σ p = 2.30) above or below the mean metacarpal and proximal phalanx radii. The model diverged for the case of R m − 0.5σ m , R p + 0.5σ p above a load of 6 N. Full size image

Figure 4b shows the parametric sweep of the joint geometry assuming constant joint separation (the initial joint eccentricity varies between 0.52 and 0.72). For a joint having the metacarpal and proximal phalanx one half standard deviation either above or below the mean11, our model yielded results comparable to the experiments, except for one case (R m − σ m , R p + σ p ). For this physiologically rare case (the radius of the articulating surfaces usually change proportionately), the model fails above an applied load of 6 N. Finally for completeness, a sweep of the other parameters including synovial fluid viscosity (which is known to change with age and disease) showed that the dynamics were relatively insensitive to these parameters as suggested by the dimensional analysis.

It is reassuring that the relatively simple model presented here yields results that closely match experimental findings, especially in the initial phase of knuckle cracking, thus underscoring the potential of numerical simulations in understanding knuckle cracking. Relaxing some of the assumptions in the model are expected to yield even better results. The model makes a number of assumptions including a 2D approximation of the joint, approximations to the pressure field, a simple acceleration profile for joint separation, an axisymmetric cavitation bubble, negligible effects of confinement on the bubble dynamics, the presence of only a single bubble, and that the acoustics are solely due to bubble collapse. One or more these assumptions, especially that of the presence of a single cavitation bubble, could very well explain why the model predictions do not closely match experiments after the dominant acoustic pressure peak. Multiple cavitation bubbles may occur in reality, but the first bubble will be expected to originate close to the centre where the lowest pressures are located. Logically, the initial dynamics will be governed by the primary bubble before interactions from other collapsing bubbles begin to influence the sound. Moreover, as the pressure drops in the centre of the joint, there will be a flow of synovial fluid from the surroundings to this low pressure region, altering the pressure field. This effect may initially be negligible but may have a significant influence on the dynamics at a later stage. The initial minor positive peak seen in the acoustic pressure measured from experiments (Fig. 3) could be from the growth phase of the bubble or from other sources6, none of which were accounted for in our model.