| Title: | Musculoskeletal Network Analysis for Physiological Data |
|---|---|
| Description: | Network science analysis of the human musculoskeletal system. Implements hypergraph representations of muscle-bone connections, physical simulation via damped harmonic oscillators, perturbation-based impact scoring, community detection, and clinical outcome prediction. Based on Murphy et al. (2018) <doi:10.1371/journal.pbio.2002811>. |
| Authors: | Yusuke Matsui [aut, cre] |
| Maintainer: | Yusuke Matsui <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.3.1 |
| Built: | 2026-10-03 13:42:26 UTC |
| Source: | https://github.com/x-biosignal/PhysioMSKNet |
CEINMS-style calibration: adjusts a small set of parameters (force scale, activation time constant, EMG non-linearity) to minimise the RMS error between the model's predicted joint moment and a measured reference moment (e.g. from inverse dynamics).
calibrateEMGDrivenModel( emg, norm_len, norm_vel, moment_arm, reference_moment, sr, params = emgDrivenParams() )calibrateEMGDrivenModel( emg, norm_len, norm_vel, moment_arm, reference_moment, sr, params = emgDrivenParams() )
emg |
EMG envelope. |
norm_len, norm_vel
|
Normalised fiber kinematics. |
moment_arm |
Moment arm (m). |
reference_moment |
Measured joint moment to match (N.m). |
sr |
Sampling rate (Hz). |
params |
Starting parameters from |
A list with calibrated params, rmse, r (correlation of predicted
vs reference moment), par (optimised values), and convergence.
sr <- 200; t <- seq(0, 0.5, by = 1 / sr) emg <- pmax(sin(2 * pi * 2 * t), 0) truth <- emgDrivenJointMoment( emgDrivenForce(emg, norm_len = 1, norm_vel = 0, sr = sr, params = emgDrivenParams(max_isometric_force = 900)), moment_arm = 0.03) fit <- calibrateEMGDrivenModel(emg, norm_len = 1, norm_vel = 0, moment_arm = 0.03, reference_moment = truth, sr = sr) fit$rmsesr <- 200; t <- seq(0, 0.5, by = 1 / sr) emg <- pmax(sin(2 * pi * 2 * t), 0) truth <- emgDrivenJointMoment( emgDrivenForce(emg, norm_len = 1, norm_vel = 0, sr = sr, params = emgDrivenParams(max_isometric_force = 900)), moment_arm = 0.03) fit <- calibrateEMGDrivenModel(emg, norm_len = 1, norm_vel = 0, moment_arm = 0.03, reference_moment = truth, sr = sr) fit$rmse
Constructs a musclePath() with a constant (angle-independent) moment arm
taken from the anatomical lookup table, i.e. a linear MTU length
L(theta) = slack_length - r * theta. Useful as a default path when no
subject-specific polynomial model is available.
defaultMusclePath(muscle, joint = NULL, slack_length = 0.3, hill = NULL)defaultMusclePath(muscle, joint = NULL, slack_length = 0.3, hill = NULL)
muscle |
Muscle name (as in |
joint |
Optional joint name (required only if the muscle appears for more than one joint). |
slack_length |
MTU slack length (m) offset (default 0.3). |
hill |
Optional Hill parameters (see |
A muscle_path object with a constant moment arm.
musclePath(), muscleTendonKinematics()
defaultMusclePath("Soleus")defaultMusclePath("Soleus")
Computes the degree distribution for bones or muscles.
degreeDistribution(hg, type = c("muscle", "bone"))degreeDistribution(hg, type = c("muscle", "bone"))
hg |
An MSKHypergraph object. |
type |
Character, either "muscle" (hyperedge degree) or "bone" (vertex degree). |
A data.frame with columns: degree, count, probability
hg <- MSKHypergraph() head(degreeDistribution(hg, type = "muscle"))hg <- MSKHypergraph() head(degreeDistribution(hg, type = "muscle"))
Detects functional communities from EMG coherence and compares them with structural communities from the MSK network.
emgCommunityCompare(pe, hg = NULL, gamma = 4.3, mapping = NULL)emgCommunityCompare(pe, hg = NULL, gamma = 4.3, mapping = NULL)
pe |
A SummarizedExperiment-like object or numeric signal matrix. |
hg |
An MSKHypergraph object (NULL loads default). |
gamma |
Resolution parameter for MSK community detection (default: 4.3). |
mapping |
Optional data.frame from |
A list with:
Named integer vector of EMG community assignments
Named integer vector of MSK community assignments
z-Rand score comparing the two partitions
The channel-to-muscle mapping used
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") emgCommunityCompare(emg, hg = hg)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") emgCommunityCompare(emg, hg = hg)
Full forward pipeline: EMG envelope -> excitation -> activation dynamics ->
Hill contractile + passive force, given normalised fiber kinematics (e.g. from
muscleTendonKinematics()).
emgDrivenForce(emg, norm_len, norm_vel, sr, params = emgDrivenParams())emgDrivenForce(emg, norm_len, norm_vel, sr, params = emgDrivenParams())
emg |
EMG envelope time series. |
norm_len |
Normalised fiber length (scalar or same length as |
norm_vel |
Normalised fiber velocity (scalar or same length as |
sr |
Sampling rate (Hz). |
params |
Parameter list from |
Muscle force time series (N).
sr <- 1000; t <- seq(0, 2, by = 1 / sr) emg <- pmax(sin(2 * pi * 1 * t), 0) * abs(rnorm(length(t), 1, 0.1)) f <- emgDrivenForce(emg, norm_len = 1, norm_vel = 0, sr = sr, params = emgDrivenParams(max_isometric_force = 800)) max(f)sr <- 1000; t <- seq(0, 2, by = 1 / sr) emg <- pmax(sin(2 * pi * 1 * t), 0) * abs(rnorm(length(t), 1, 0.1)) f <- emgDrivenForce(emg, norm_len = 1, norm_vel = 0, sr = sr, params = emgDrivenParams(max_isometric_force = 800)) max(f)
Forward pipeline as emgDrivenForce() but with tendon compliance: at each
time the muscle-tendon force equilibrium is solved for the fiber length
(quasi-static), so tendon stretch is accounted for. Requires the MTU length
trajectory (e.g. from muscleTendonKinematics()).
emgDrivenForceElastic( emg, mtu_length, sr, optimal_fiber_length, tendon_slack_length, params = emgDrivenParams() )emgDrivenForceElastic( emg, mtu_length, sr, optimal_fiber_length, tendon_slack_length, params = emgDrivenParams() )
emg |
EMG envelope. |
mtu_length |
MTU length (m); scalar or same length as |
sr |
Sampling rate (Hz). |
optimal_fiber_length, tendon_slack_length
|
Muscle geometry (m). |
params |
Parameter list from |
Muscle-tendon force time series (N).
sr <- 500; t <- seq(0, 1, by = 1 / sr) emg <- pmax(sin(2 * pi * 2 * t), 0) f <- emgDrivenForceElastic(emg, mtu_length = 0.30, sr = sr, optimal_fiber_length = 0.10, tendon_slack_length = 0.19) max(f)sr <- 500; t <- seq(0, 1, by = 1 / sr) emg <- pmax(sin(2 * pi * 2 * t), 0) f <- emgDrivenForceElastic(emg, mtu_length = 0.30, sr = sr, optimal_fiber_length = 0.10, tendon_slack_length = 0.19) max(f)
Combines muscle forces with moment arms into a net joint moment (N.m).
emgDrivenJointMoment(forces, moment_arm)emgDrivenJointMoment(forces, moment_arm)
forces |
A force time series (one muscle) or a time x muscle matrix. |
moment_arm |
Moment arm(s) in m: a scalar/vector matching |
Net joint moment time series (N.m).
emgDrivenJointMoment(c(100, 150, 200), moment_arm = 0.03)emgDrivenJointMoment(c(100, 150, 200), moment_arm = 0.03)
Default EMG-driven muscle parameters
emgDrivenParams( max_isometric_force = 1000, pennation = 0, tau_act = 0.01, tau_deact = 0.04, emg_nonlin = 0, max_emg = NULL )emgDrivenParams( max_isometric_force = 1000, pennation = 0, tau_act = 0.01, tau_deact = 0.04, emg_nonlin = 0, max_emg = NULL )
max_isometric_force |
Peak isometric force |
pennation |
Pennation angle in radians (default 0). |
tau_act, tau_deact
|
Activation dynamics time constants (s). |
emg_nonlin |
EMG-to-activation non-linearity coefficient. |
max_emg |
Optional EMG normaliser. |
A named list of parameters.
emgDrivenParams(max_isometric_force = 800)emgDrivenParams(max_isometric_force = 800)
Tests whether EMG activation levels differ across MSK community assignments using a Kruskal-Wallis test.
emgMSKEnrichment(pe, hg = NULL, gamma = 4.3, mapping = NULL)emgMSKEnrichment(pe, hg = NULL, gamma = 4.3, mapping = NULL)
pe |
A SummarizedExperiment-like object or numeric signal matrix. |
hg |
An MSKHypergraph object (NULL loads default). |
gamma |
Resolution parameter for MSK community detection (default: 4.3). |
mapping |
Optional data.frame from |
A list with:
Data frame with community, mean/median activation
Kruskal-Wallis test result
Named numeric vector of RMS activation
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") emgMSKEnrichment(emg, hg = hg)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") emgMSKEnrichment(emg, hg = hg)
Computes EMG functional coherence and compares it with structural adjacency from the MSK muscle graph using a Mantel test.
emgStructuralCoherence(pe, hg = NULL, freq_band = c(20, 50), mapping = NULL)emgStructuralCoherence(pe, hg = NULL, freq_band = c(20, 50), mapping = NULL)
pe |
A SummarizedExperiment-like object or a numeric signal matrix (time x channels). |
hg |
An MSKHypergraph object (NULL loads default). |
freq_band |
Numeric vector of length 2, frequency band in Hz for coherence (default: c(20, 50) for EMG beta/gamma). |
mapping |
Optional data.frame from |
A list with:
EMG functional coherence matrix
MSK structural adjacency (matched subset)
Mantel correlation coefficient
Permutation-based p-value
Names of matched muscles
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") emgStructuralCoherence(emg, hg = hg, freq_band = c(20, 50))C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") emgStructuralCoherence(emg, hg = hg, freq_band = c(20, 50))
Rectifies and normalises an EMG envelope to [0, 1], with an optional
A-shape (exponential) EMG-to-activation non-linearity.
emgToExcitation(emg, max_emg = NULL, nonlin = 0)emgToExcitation(emg, max_emg = NULL, nonlin = 0)
emg |
EMG envelope (already rectified/smoothed, or raw — abs is taken). |
max_emg |
Normalising maximum (default: the signal maximum). |
nonlin |
Non-linearity coefficient |
Excitation in [0, 1].
emgToExcitation(c(0.1, 0.5, 0.9, 0.3))emgToExcitation(c(0.1, 0.5, 0.9, 0.3))
Matches EMG channel names from a PhysioExperiment/SummarizedExperiment object to muscles in an MSK hypergraph.
emgToMSKMapping(pe, hg = NULL, method = c("exact", "fuzzy"), threshold = 0.8)emgToMSKMapping(pe, hg = NULL, method = c("exact", "fuzzy"), threshold = 0.8)
pe |
A SummarizedExperiment-like object with channel names in
|
hg |
An MSKHypergraph object (NULL loads default). |
method |
Character, matching method: "exact" or "fuzzy". |
threshold |
Numeric, fuzzy matching threshold (default: 0.8). |
A data.frame with columns: channel_idx, channel_name, muscle_idx, muscle_name, match_quality.
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") mapping <- emgToMSKMapping(emg, method = "fuzzy") mappingset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") mapping <- emgToMSKMapping(emg, method = "fuzzy") mapping
Solves the quasi-static (isometric) muscle-tendon force balance for the fiber
length: the fiber active+passive force along the tendon equals the tendon
force, with tendon length = MTU length - fiber length * cos(pennation).
equilibriumFiberLength( activation, mtu_length, optimal_fiber_length, tendon_slack_length, max_isometric_force = 1000, pennation = 0, e0t = 0.04 )equilibriumFiberLength( activation, mtu_length, optimal_fiber_length, tendon_slack_length, max_isometric_force = 1000, pennation = 0, e0t = 0.04 )
activation |
Activation in |
mtu_length |
Muscle-tendon-unit length (m). |
optimal_fiber_length |
Optimal fiber length (m). |
tendon_slack_length |
Tendon slack length (m). |
max_isometric_force |
Peak isometric force (N). |
pennation |
Pennation angle (rad). |
e0t |
Tendon reference strain (default 0.04). |
A list with fiber_length, tendon_length, norm_fiber_length,
norm_tendon_length, and force (N).
equilibriumFiberLength(0.5, mtu_length = 0.30, optimal_fiber_length = 0.10, tendon_slack_length = 0.20)equilibriumFiberLength(0.5, mtu_length = 0.30, optimal_fiber_length = 0.10, tendon_slack_length = 0.20)
First-order activation dynamics (Thelen 2003): activation lags excitation, with faster activation than deactivation.
excitationToActivation( excitation, dt, tau_act = 0.01, tau_deact = 0.04, a0 = NULL )excitationToActivation( excitation, dt, tau_act = 0.01, tau_deact = 0.04, a0 = NULL )
excitation |
Neural excitation in |
dt |
Time step in seconds. |
tau_act |
Activation time constant (default 0.010 s). |
tau_deact |
Deactivation time constant (default 0.040 s). |
a0 |
Initial activation (default: |
Activation time series in [0, 1].
u <- c(0, 0, 1, 1, 1, 0, 0) excitationToActivation(u, dt = 0.01)u <- c(0, 0, 1, 1, 1, 0, 0) excitationToActivation(u, dt = 0.01)
Gaussian active force-length curve, peaking at the optimal fiber length (normalised length 1).
forceLengthActive(norm_len, width = 0.45)forceLengthActive(norm_len, width = 0.45)
norm_len |
Normalised fiber length (optimal = 1). |
width |
Gaussian width parameter (default 0.45). |
Multiplier in [0, 1].
forceLengthActive(seq(0.6, 1.4, by = 0.2))forceLengthActive(seq(0.6, 1.4, by = 0.2))
Exponential passive element (Thelen 2003): zero below optimal length, rising as the fiber is stretched beyond it.
forceLengthPassive(norm_len, k_pe = 4, e0 = 0.6)forceLengthPassive(norm_len, k_pe = 4, e0 = 0.6)
norm_len |
Normalised fiber length. |
k_pe |
Exponential shape (default 4). |
e0 |
Passive strain at maximum isometric force (default 0.6). |
Passive multiplier (>= 0).
forceLengthPassive(c(1.0, 1.1, 1.2))forceLengthPassive(c(1.0, 1.1, 1.2))
Normalised Hill force-velocity relation. Velocity is normalised fiber velocity with negative = shortening (concentric), positive = lengthening (eccentric).
forceVelocity(norm_vel, af = 0.25, f_ecc = 1.8, k_ecc = 0.15)forceVelocity(norm_vel, af = 0.25, f_ecc = 1.8, k_ecc = 0.15)
norm_vel |
Normalised fiber velocity (concentric < 0). |
af |
Hill shape parameter for shortening (default 0.25). |
f_ecc |
Eccentric force asymptote (default 1.8). |
k_ecc |
Eccentric saturation constant (default 0.15). |
Multiplier (0 at max shortening, 1 at isometric, up to f_ecc).
forceVelocity(c(-0.5, 0, 0.5))forceVelocity(c(-0.5, 0, 0.5))
Total muscle-tendon force from activation and normalised fiber kinematics:
.
hillMuscleForce( activation, norm_len, norm_vel, max_isometric_force = 1000, pennation = 0 )hillMuscleForce( activation, norm_len, norm_vel, max_isometric_force = 1000, pennation = 0 )
activation |
Activation in |
norm_len |
Normalised fiber length. |
norm_vel |
Normalised fiber velocity. |
max_isometric_force |
Peak isometric force (N). |
pennation |
Pennation angle (rad). |
Muscle force in N.
hillMuscleForce(activation = 0.5, norm_len = 1, norm_vel = 0, max_isometric_force = 800)hillMuscleForce(activation = 0.5, norm_len = 1, norm_vel = 0, max_isometric_force = 800)
Returns the degree of each hyperedge (muscle), i.e., the number of bones each muscle attaches to.
hyperedgeDegree(hg)hyperedgeDegree(hg)
hg |
An MSKHypergraph object or incidence matrix. |
Named numeric vector of muscle degrees.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hyperedgeDegree(MSKHypergraph(C))C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hyperedgeDegree(MSKHypergraph(C))
Analyzes movement synchrony from IMU sensors in the context of MSK network communities. Computes within-community vs between-community sensor synchrony and tests for significance.
imuCommunityDynamics( imu_data, hg = NULL, gamma = 4.3, mapping = NULL, window_sec = NULL )imuCommunityDynamics( imu_data, hg = NULL, gamma = 4.3, mapping = NULL, window_sec = NULL )
imu_data |
Named list of per-sensor data (see |
hg |
An MSKHypergraph object (NULL loads default). |
gamma |
Numeric, resolution parameter for community detection (default 4.3). |
mapping |
Pre-computed mapping (optional). |
window_sec |
Numeric or NULL. If provided, computes time-resolved
synchrony in sliding windows of this duration (seconds). Requires
|
A list with:
Mean synchrony within communities
Mean synchrony between communities
Within/between ratio
Wilcoxon test p-value
Data frame if window_sec provided, NULL otherwise
set.seed(1) n <- 200 imu <- list( thigh = matrix(rnorm(n * 3), n, 3), shank = matrix(rnorm(n * 3), n, 3), foot = matrix(rnorm(n * 3), n, 3), pelvis = matrix(rnorm(n * 3), n, 3)) res <- imuCommunityDynamics(imu) res$ratioset.seed(1) n <- 200 imu <- list( thigh = matrix(rnorm(n * 3), n, 3), shank = matrix(rnorm(n * 3), n, 3), foot = matrix(rnorm(n * 3), n, 3), pelvis = matrix(rnorm(n * 3), n, 3)) res <- imuCommunityDynamics(imu) res$ratio
Computes bone stress from IMU acceleration/angular velocity data and propagates it through the MSK incidence matrix to predict muscle vulnerability.
imuImpactPrediction( imu_data, hg = NULL, mapping = NULL, stress_metric = c("acceleration", "angular_velocity", "jerk", "composite"), use_proxy = TRUE )imuImpactPrediction( imu_data, hg = NULL, mapping = NULL, stress_metric = c("acceleration", "angular_velocity", "jerk", "composite"), use_proxy = TRUE )
imu_data |
A named list of per-sensor data, or a matrix (see
|
hg |
An MSKHypergraph object (NULL loads default). |
mapping |
Pre-computed mapping (optional). |
stress_metric |
Character: |
use_proxy |
Logical, use degree-based proxy for impact deviation (faster) or run full simulation. |
A list with:
Named numeric vector of muscle vulnerability scores
Named numeric vector of bone stress
Named numeric of stress propagated to muscles
Data frame ranked by vulnerability (descending)
set.seed(1) n <- 200 imu <- list( thigh = matrix(rnorm(n * 3), n, 3), shank = matrix(rnorm(n * 3), n, 3), foot = matrix(rnorm(n * 3), n, 3)) pred <- imuImpactPrediction(imu, stress_metric = "acceleration") head(pred$ranking)set.seed(1) n <- 200 imu <- list( thigh = matrix(rnorm(n * 3), n, 3), shank = matrix(rnorm(n * 3), n, 3), foot = matrix(rnorm(n * 3), n, 3)) pred <- imuImpactPrediction(imu, stress_metric = "acceleration") head(pred$ranking)
Computes pairwise kinematic coupling from IMU orientation or acceleration data and compares it to the structural adjacency of the MSK bone graph using a Mantel test.
imuNetworkKinematics( imu_data, hg = NULL, mapping = NULL, signal = c("orientation", "acceleration"), method = c("correlation", "mutual_info") )imuNetworkKinematics( imu_data, hg = NULL, mapping = NULL, signal = c("orientation", "acceleration"), method = c("correlation", "mutual_info") )
imu_data |
A named list of per-sensor data. Each element should be a
matrix or data.frame with orientation/acceleration columns. The list names
are sensor placement names (e.g., |
hg |
An MSKHypergraph object (NULL loads default). |
mapping |
A pre-computed mapping from |
signal |
Character, which signal to use for coupling:
|
method |
Character, coupling method: |
A list with:
Pairwise coupling matrix between mapped bones
Corresponding MSK bone graph adjacency
Mantel test correlation coefficient
Permutation p-value
Data frame of sensor-to-bone mapping used
set.seed(1) n <- 200 imu <- list( thigh = matrix(rnorm(n * 3), n, 3), shank = matrix(rnorm(n * 3), n, 3), foot = matrix(rnorm(n * 3), n, 3), pelvis = matrix(rnorm(n * 3), n, 3)) res <- imuNetworkKinematics(imu, signal = "orientation") res$mapped_sensorsset.seed(1) n <- 200 imu <- list( thigh = matrix(rnorm(n * 3), n, 3), shank = matrix(rnorm(n * 3), n, 3), foot = matrix(rnorm(n * 3), n, 3), pelvis = matrix(rnorm(n * 3), n, 3)) res <- imuNetworkKinematics(imu, signal = "orientation") res$mapped_sensors
Matches IMU sensor placement names to bones in the MSK hypergraph using a curated lookup table with fuzzy matching fallback. Supports naming conventions from Xsens, APDM, Shimmer, and generic body-segment labels.
imuToMSKMapping( sensors, hg = NULL, method = c("exact", "fuzzy"), threshold = 0.7 )imuToMSKMapping( sensors, hg = NULL, method = c("exact", "fuzzy"), threshold = 0.7 )
sensors |
Character vector of IMU sensor placement names
(e.g., |
hg |
An MSKHypergraph object (NULL loads default). |
method |
Character, matching method: "exact" or "fuzzy". |
threshold |
Numeric, fuzzy matching threshold (default: 0.7). |
A data.frame with columns: sensor_name, bone_idx, bone_name, match_quality, match_method.
mapping <- imuToMSKMapping(c("upper_arm", "thigh", "lumbar")) mappingmapping <- imuToMSKMapping(c("upper_arm", "thigh", "lumbar")) mapping
Loads the bipartite incidence matrix C where rows are bones (vertices) and columns are muscles (hyperedges). Entry C[i,j] = 1 indicates that muscle j attaches to bone i.
loadIncidenceMatrix()loadIncidenceMatrix()
A sparse Matrix (dgCMatrix) of dimensions 173 x 270
Murphy AC et al. (2018) PLOS Biology.
C <- loadIncidenceMatrix() dim(C)C <- loadIncidenceMatrix() dim(C)
Loads the complete musculoskeletal network dataset from Murphy et al. (2018). Returns a list containing the incidence matrix, muscle metadata, and validation data for reproducing the paper's results.
loadMSKData()loadMSKData()
A list with components:
Sparse incidence matrix C (bones x muscles)
Data frame with muscle names, community assignments, homunculus categories
Character vector of bone/vertex names
Character vector of muscle/hyperedge names
Murphy AC et al. (2018) "Structure, function, and control of the human musculoskeletal network." PLOS Biology 16(1): e2002811.
data <- loadMSKData() dim(data$incidence) head(data$bone_names)data <- loadMSKData() dim(data$incidence) head(data$bone_names)
Loads metadata for all 270 muscles including community assignments and homunculus category labels from the paper.
loadMuscleMetadata()loadMuscleMetadata()
A data.frame with columns: index, muscle, community, homunculus_category
Murphy AC et al. (2018) PLOS Biology.
meta <- loadMuscleMetadata() head(meta)meta <- loadMuscleMetadata() head(meta)
Loads validation datasets used to reproduce key figures from the paper.
loadValidationData( dataset = c("degree_distribution", "impact_vs_recovery", "homunculus_deviation", "fmri_activation", "homunculus_coordinates", "impact_vs_path", "impact_scores") )loadValidationData( dataset = c("degree_distribution", "impact_vs_recovery", "homunculus_deviation", "fmri_activation", "homunculus_coordinates", "impact_vs_path", "impact_scores") )
dataset |
Character string specifying which dataset to load:
|
A data.frame with the requested validation data
Murphy AC et al. (2018) PLOS Biology.
rec <- loadValidationData("impact_vs_recovery") head(rec)rec <- loadValidationData("impact_vs_recovery") head(rec)
Computes displacement synchrony (correlation) between matched bones and tests whether within-community pairs are more synchronized than between-community pairs using a Wilcoxon test.
mocapCommunityDynamics( pe_mocap, hg = NULL, gamma = 4.3, mapping = NULL, window_sec = NULL )mocapCommunityDynamics( pe_mocap, hg = NULL, gamma = 4.3, mapping = NULL, window_sec = NULL )
pe_mocap |
A numeric matrix (frames x segments) with MoCap data. |
hg |
An MSKHypergraph object (NULL loads default). |
gamma |
Resolution parameter for community detection (default: 4.3). |
mapping |
Optional data.frame from |
window_sec |
Optional numeric, window size in seconds for time-resolved analysis (NULL for global analysis only). |
A list with:
Mean within-community synchrony
Mean between-community synchrony
Within/between ratio
Wilcoxon test p-value
Data frame with per-window results (if window_sec set)
set.seed(6) seg <- matrix(cumsum(rnorm(300 * 3)), 300, 3) colnames(seg) <- c("humerus", "radius", "femur") attr(seg, "sr") <- 120 mocapCommunityDynamics(seg)set.seed(6) seg <- matrix(cumsum(rnorm(300 * 3)), 300, 3) colnames(seg) <- c("humerus", "radius", "femur") attr(seg, "sr") <- 120 mocapCommunityDynamics(seg)
Computes kinematic stress per bone from MoCap data and propagates it through the MSK incidence matrix to estimate muscle vulnerability.
mocapImpactPrediction( pe_mocap, hg = NULL, mapping = NULL, stress_metric = c("acceleration", "jerk", "range"), use_proxy = TRUE )mocapImpactPrediction( pe_mocap, hg = NULL, mapping = NULL, stress_metric = c("acceleration", "jerk", "range"), use_proxy = TRUE )
pe_mocap |
A numeric matrix (frames x segments) with MoCap data. |
hg |
An MSKHypergraph object (NULL loads default). |
mapping |
Optional data.frame from |
stress_metric |
Character, kinematic stress metric: "acceleration", "jerk", or "range". |
use_proxy |
Logical, if TRUE uses degree-based proxy instead of full simulation for impact deviation (default: TRUE). |
A list with:
Named numeric vector of muscle vulnerability scores
Named numeric vector of bone stress values
Named numeric vector of muscle stress exposure
Data frame ranking muscles by vulnerability
set.seed(6) seg <- matrix(cumsum(rnorm(300 * 3)), 300, 3) colnames(seg) <- c("humerus", "radius", "femur") attr(seg, "sr") <- 120 mocapImpactPrediction(seg)set.seed(6) seg <- matrix(cumsum(rnorm(300 * 3)), 300, 3) colnames(seg) <- c("humerus", "radius", "femur") attr(seg, "sr") <- 120 mocapImpactPrediction(seg)
Computes pairwise movement coupling between MoCap segments and compares the kinematic coupling matrix with MSK structural adjacency.
mocapNetworkKinematics( pe_mocap, hg = NULL, mapping = NULL, method = c("correlation", "mutual_info") )mocapNetworkKinematics( pe_mocap, hg = NULL, mapping = NULL, method = c("correlation", "mutual_info") )
pe_mocap |
A numeric matrix (frames x segments) or SummarizedExperiment with MoCap position data. |
hg |
An MSKHypergraph object (NULL loads default). |
mapping |
Optional data.frame from |
method |
Character, coupling method: "correlation" or "mutual_info". |
A list with:
Pairwise coupling matrix
MSK bone adjacency (matched subset)
Mantel correlation coefficient
Permutation-based p-value
set.seed(6) seg <- matrix(cumsum(rnorm(300 * 3)), 300, 3) colnames(seg) <- c("humerus", "radius", "femur") attr(seg, "sr") <- 120 mocapNetworkKinematics(seg)set.seed(6) seg <- matrix(cumsum(rnorm(300 * 3)), 300, 3) colnames(seg) <- c("humerus", "radius", "femur") attr(seg, "sr") <- 120 mocapNetworkKinematics(seg)
Matches MoCap segment names from a SkeletonModel or character vector to bones in an MSK hypergraph using a curated lookup table with fuzzy matching fallback.
mocapToMSKMapping( skeleton, hg = NULL, method = c("exact", "fuzzy"), threshold = 0.7 )mocapToMSKMapping( skeleton, hg = NULL, method = c("exact", "fuzzy"), threshold = 0.7 )
skeleton |
A PhysioMoCap SkeletonModel object or character vector of segment names. |
hg |
An MSKHypergraph object (NULL loads default). |
method |
Character, matching method: "exact" or "fuzzy". |
threshold |
Numeric, fuzzy matching threshold (default: 0.7). |
A data.frame with columns: segment_name, bone_idx, bone_name, match_quality, match_method.
mapping <- mocapToMSKMapping(c("humerus", "radius", "femur")) mappingmapping <- mocapToMSKMapping(c("humerus", "radius", "femur")) mapping
Returns muscle moment arms. With no model, this is the anatomical lookup
table promoted from the package's internal reference values (signed by the
agonist/antagonist direction); with a musclePath() model, the moment arm
is the joint-angle-dependent value -dL/dtheta evaluated at angles.
momentArm(muscle = NULL, joint = NULL, model = NULL, angles = NULL)momentArm(muscle = NULL, joint = NULL, model = NULL, angles = NULL)
muscle |
Optional muscle name to filter the lookup table (case
insensitive). Ignored when |
joint |
Optional joint name to filter the lookup table. |
model |
Optional |
angles |
Named numeric vector of joint angles in radians (one per joint
the |
Without model: a data frame with muscle_name, joint_name,
moment_arm_m (magnitude), direction, and signed_moment_arm_m. With
model: a named numeric vector of signed moment arms (m) per joint.
Sherman MA, Seth A, Delp SL (2013). "What is a moment arm? Calculating muscle effectiveness in biomechanical models using generalized coordinates."
musclePath(), muscleTendonKinematics()
momentArm("Gastrocnemius") mp <- musclePath("ankle", coefficients = list(ankle = c(-0.05))) momentArm(model = mp, angles = c(ankle = 0.2))momentArm("Gastrocnemius") mp <- musclePath("ankle", coefficients = list(ankle = c(-0.05))) momentArm(model = mp, angles = c(ankle = 0.2))
Modifies a rehabilitation protocol based on patient reassessment results. Uses decision rules to determine whether to advance, continue, reduce, or modify the current protocol.
mskAdaptProtocol(reassessment, current_protocol, hg = NULL)mskAdaptProtocol(reassessment, current_protocol, hg = NULL)
reassessment |
An |
current_protocol |
An |
hg |
An MSKHypergraph object (NULL loads default). |
An S3 object of class "MSKAdaptedProtocol" with:
character: "advance", "continue", "reduce", "modify"
data.frame with modified exercise prescription
character explanation
character vector of muscles needing attention
exercises to drop
new exercises to add
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") reassessment <- mskReassess(list(emg = emg), outcome, hg = hg) protocol <- mskRehabProtocol("Biceps Brachii", hg = hg) adapted <- mskAdaptProtocol(reassessment, protocol, hg = hg) print(adapted)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") reassessment <- mskReassess(list(emg = emg), outcome, hg = hg) protocol <- mskRehabProtocol("Biceps Brachii", hg = hg) adapted <- mskAdaptProtocol(reassessment, protocol, hg = hg) print(adapted)
Attaches anatomical annotations from PhysioAnnotationHub to an MSK hypergraph, matching muscle and bone names between the two data sources. Annotations include body region, innervation, actions, and spinal levels.
mskAnnotate(hg = NULL, hub = NULL)mskAnnotate(hg = NULL, hub = NULL)
hg |
An MSKHypergraph object. If NULL, loads default 173-bone/270-muscle network. |
hub |
A PhysioAnnotationHub object. If NULL, loads via
|
An annotated MSKHypergraph with additional fields:
Data frame of muscle annotations merged from hub
Data frame of bone annotations merged from hub
Logical flag indicating annotations are attached
List with muscle and bone match rates
Annotations are derived from PhysioAnnotationHub's curated knowledge graph. Name matching uses fuzzy matching with a 0.2 edit distance threshold, which may produce incorrect matches for similarly named structures. Always verify critical annotations against primary anatomical references.
Murphy AC et al. (2018) PLOS Biology 16(1): e2002811.
if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { hg <- mskAnnotate() head(hg$muscle_annotations) }if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { hg <- mskAnnotate() head(hg$muscle_annotations) }
Computes betweenness centrality for the projected graph. Uses igraph if available, otherwise a pure-R BFS implementation.
mskBetweenness(hg, type = c("bone", "muscle"))mskBetweenness(hg, type = c("bone", "muscle"))
hg |
An MSKHypergraph object. |
type |
Character, "bone" or "muscle" projection. |
Named numeric vector of betweenness centrality values.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) mskBetweenness(MSKHypergraph(C), type = "bone")C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) mskBetweenness(MSKHypergraph(C), type = "bone")
Generates a clinical evidence report for a set of injured muscles by querying the knowledge graph for ICD-10 codes, ICF codes, innervation, functional impact, and related muscles (synergists and antagonists).
mskClinicalEvidence(injury_muscles, hub = NULL, hg = NULL)mskClinicalEvidence(injury_muscles, hub = NULL, hg = NULL)
injury_muscles |
Character or integer vector identifying injured muscles. |
hub |
A PhysioAnnotationHub object (NULL loads default). |
hg |
An MSKHypergraph object (NULL loads default). |
An S3 object of class "MSKClinicalEvidence" with:
Data frame of injury muscle annotations
Matching ICD-10 entries
Matching ICF entries
Nerves innervating injured muscles
Spinal segments involved
Affected actions/movements
Synergistic muscles from KG
Antagonistic muscles from KG
ICD-10 and ICF code mappings are based on knowledge graph associations and may not capture all valid codes for a clinical scenario. This tool provides evidence aggregation for research; clinical coding should follow local guidelines and be reviewed by qualified professionals.
if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { evidence <- mskClinicalEvidence(c("Biceps Brachii", "Deltoid")) print(evidence) }if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { evidence <- mskClinicalEvidence(c("Biceps Brachii", "Deltoid")) print(evidence) }
Predicts recovery time, identifies compensatory muscles, and estimates secondary injury risk based on MSK network impact analysis.
mskClinicalPredictor(injury_muscles, hg = NULL, sim = NULL, verbose = TRUE)mskClinicalPredictor(injury_muscles, hg = NULL, sim = NULL, verbose = TRUE)
injury_muscles |
Character or integer vector identifying injured muscles. |
hg |
An MSKHypergraph object (NULL loads default 173-bone/270-muscle network). |
sim |
An MSKSimulation object (NULL creates one with default parameters). |
verbose |
Logical, print progress messages (default: TRUE). |
An S3 object of class "MSKClinicalPrediction" with:
Data frame with predicted recovery weeks and CI per muscle
List of compensatory muscles per injured muscle
Data frame of secondary injury risk scores
Resolved muscle names
Resolved integer indices
The recovery model was validated on 14 aggregate muscle groups, not individual muscles. Patient factor adjustments are heuristic, not independently validated. This is a research exploration tool, not a clinical diagnostic.
Murphy AC et al. (2018) PLOS Biology 16(1): e2002811.
# Build a small musculoskeletal hypergraph (bones x muscles incidence) C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) pred <- mskClinicalPredictor("Biceps Brachii", hg = hg, verbose = FALSE) print(pred)# Build a small musculoskeletal hypergraph (bones x muscles incidence) C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) pred <- mskClinicalPredictor("Biceps Brachii", hg = hg, verbose = FALSE) print(pred)
Closeness Centrality
mskCloseness(hg, type = c("bone", "muscle"))mskCloseness(hg, type = c("bone", "muscle"))
hg |
An MSKHypergraph object. |
type |
Character, "bone" or "muscle" projection. |
Named numeric vector of closeness centrality values.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) mskCloseness(MSKHypergraph(C), type = "bone")C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) mskCloseness(MSKHypergraph(C), type = "bone")
Detects communities in the projected muscle-muscle graph using the Louvain algorithm with a resolution parameter gamma.
mskCommunityDetect(hg = NULL, gamma = 4.3, type = c("muscle", "bone"))mskCommunityDetect(hg = NULL, gamma = 4.3, type = c("muscle", "bone"))
hg |
An MSKHypergraph object. If NULL, loads built-in data. |
gamma |
Resolution parameter for modularity (default: 4.3, as in paper). Higher values produce more, smaller communities. |
type |
Projection type: "muscle" (default, as in paper) or "bone". |
A list with:
Named integer vector of community assignments
Number of communities detected
Modularity value Q
Resolution parameter used
Table of community sizes
Murphy AC et al. (2018) PLOS Biology 16(1): e2002811.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) cm <- mskCommunityDetect(MSKHypergraph(C), gamma = 1) cm$membershipC <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) cm <- mskCommunityDetect(MSKHypergraph(C), gamma = 1) cm$membership
Generates a detailed functional profile of a musculoskeletal community by combining network topology metrics with knowledge graph annotations. This enables characterization of communities by their anatomical region, innervation patterns, primary actions, and spinal segment involvement.
mskCommunityProfile(hg = NULL, community_id, hub = NULL, gamma = 4.3)mskCommunityProfile(hg = NULL, community_id, hub = NULL, gamma = 4.3)
hg |
An MSKHypergraph object. If NULL, loads default network. |
community_id |
Integer, the community to profile (1-based). |
hub |
A PhysioAnnotationHub object (NULL loads default). |
gamma |
Numeric, resolution parameter for community detection (default: 4.3, as in Murphy et al. 2018). |
An S3 object of class "MSKCommunityProfile" with:
Integer, the profiled community
Character vector of muscle names in the community
Integer, number of muscles
Table of primary action distribution
Table of innervating nerves distribution
Table of body region distribution
Table of spinal level distribution
Most frequent primary action
Most frequent innervating nerve
Most frequent body region
Mean hyperedge degree of community muscles
Mean impact deviation score
Community membership depends on the resolution parameter gamma and the stochastic Louvain algorithm. Profile annotations depend on KG completeness. Community boundaries are network-derived, not anatomical boundaries.
Murphy AC et al. (2018) PLOS Biology 16(1): e2002811.
if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) profile <- mskCommunityProfile(community_id = 1, hg = hg) print(profile) }if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) profile <- mskCommunityProfile(community_id = 1, hg = hg) print(profile) }
Tracks changes in compensation patterns across multiple timepoints, identifying onset, resolution, and trends.
mskCompensationEvolution( timepoints_emg, injured_muscles, hg = NULL, emg_mapping = NULL, sr = NULL, z_threshold = 1.96 )mskCompensationEvolution( timepoints_emg, injured_muscles, hg = NULL, emg_mapping = NULL, sr = NULL, z_threshold = 1.96 )
timepoints_emg |
Named list of EMG matrices. The first element is treated as baseline. |
injured_muscles |
Character vector of injured muscle names or integer indices. |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
sr |
Optional sampling rate override. |
z_threshold |
Numeric, z-score threshold (default: 1.96). |
An S3 object of class "MSKCompensationEvolution" with:
Data.frame of per-timepoint per-muscle z-scores
Per-muscle first timepoint of compensation
Per-muscle first timepoint of resolution
Per-muscle trend classification
Data.frame of per-timepoint summary statistics
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tps <- list(week0 = mk(1), week2 = mk(2), week4 = mk(3)) evolution <- mskCompensationEvolution(tps, "Biceps Brachii", hg = hg) evolutionC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tps <- list(week0 = mk(1), week2 = mk(2), week4 = mk(3)) evolution <- mskCompensationEvolution(tps, "Biceps Brachii", hg = hg) evolution
Constructs a compensation-weighted subgraph from compensation detection results, identifying compensation chains and hub muscles.
mskCompensationNetwork(compensation_result, hg = NULL, emg_mapping = NULL)mskCompensationNetwork(compensation_result, hg = NULL, emg_mapping = NULL)
compensation_result |
An |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
A list with:
Weighted compensation adjacency matrix
List of compensation chains
Muscles in multiple compensation chains
MSK communities affected by compensation
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- mskDetectCompensation(emg, emg_base, "Biceps Brachii", hg = hg) net <- mskCompensationNetwork(result, hg = hg) netC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- mskDetectCompensation(emg, emg_base, "Biceps Brachii", hg = hg) net <- mskCompensationNetwork(result, hg = hg) net
Computes risk scores for compensating muscles based on biomechanical overuse principles, integrating network topology, duration of compensation, and loading intensity.
mskCompensationRiskScore( compensation_result, hg = NULL, duration_weeks = 0, load_intensity = c("low", "moderate", "high") )mskCompensationRiskScore( compensation_result, hg = NULL, duration_weeks = 0, load_intensity = c("low", "moderate", "high") )
compensation_result |
An |
hg |
An MSKHypergraph object (NULL loads default). |
duration_weeks |
Numeric, how long compensation has been occurring. |
load_intensity |
Character, one of "low", "moderate", "high". |
An S3 object of class "MSKCompensationRisk" with:
Data.frame of per-muscle risk scores
Numeric overall risk score
Character risk category
Name of highest-risk muscle
Clinical action recommendation
Per compensating muscle: risk = z_excess * degree_normalized * duration_factor * load_factor where z_excess = max(0, z_score - z_threshold), degree_normalized = hyperedgeDegree / mean_degree, duration_factor = 1 + log(1 + duration_weeks), load_factor = intensity multiplier (0.5, 1.0, 1.5).
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- mskDetectCompensation(emg, emg_base, "Biceps Brachii", hg = hg) risk <- mskCompensationRiskScore(result, hg = hg, duration_weeks = 4, load_intensity = "moderate") riskC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- mskDetectCompensation(emg, emg_base, "Biceps Brachii", hg = hg) risk <- mskCompensationRiskScore(result, hg = hg, duration_weeks = 4, load_intensity = "moderate") risk
Orchestrator that runs all compensation analyses and returns a unified summary. Uses tryCatch for each sub-analysis so partial results are available even if some analyses fail.
mskCompensationSummary( emg, emg_baseline, injured_muscles, timepoints_emg = NULL, hg = NULL, sr = NULL, z_threshold = 1.96, duration_weeks = 0, load_intensity = c("low", "moderate", "high") )mskCompensationSummary( emg, emg_baseline, injured_muscles, timepoints_emg = NULL, hg = NULL, sr = NULL, z_threshold = 1.96, duration_weeks = 0, load_intensity = c("low", "moderate", "high") )
emg |
Current EMG data (matrix or SummarizedExperiment). |
emg_baseline |
Baseline/pre-injury EMG data. |
injured_muscles |
Character vector of injured muscle names or integer indices. |
timepoints_emg |
Optional named list of EMG matrices for evolution analysis. |
hg |
An MSKHypergraph object (NULL loads default). |
sr |
Optional sampling rate override. |
z_threshold |
Numeric, z-score threshold (default: 1.96). |
duration_weeks |
Numeric, compensation duration for risk scoring. |
load_intensity |
Character, load intensity for risk scoring. |
An S3 object of class "MSKCompensationSummary" with
sub-results and available_analyses vector.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") s <- mskCompensationSummary(emg, emg_base, "Biceps Brachii", hg = hg, duration_weeks = 4, load_intensity = "moderate") sC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") s <- mskCompensationSummary(emg, emg_base, "Biceps Brachii", hg = hg, duration_weeks = 4, load_intensity = "moderate") s
Runs Louvain community detection multiple times and extracts a consensus partition. The paper uses 100 runs and takes the most frequent partition for each node pair.
mskConsensusPartition( hg = NULL, gamma = 4.3, n_runs = 100L, type = c("muscle", "bone") )mskConsensusPartition( hg = NULL, gamma = 4.3, n_runs = 100L, type = c("muscle", "bone") )
hg |
An MSKHypergraph object. |
gamma |
Resolution parameter (default: 4.3). |
n_runs |
Number of runs for consensus (default: 100). |
type |
Projection type: "muscle" or "bone". |
Same structure as mskCommunityDetect, but with consensus partition.
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) cp <- mskConsensusPartition(MSKHypergraph(C), gamma = 1, n_runs = 5) cp$n_communitiesC <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) cp <- mskConsensusPartition(MSKHypergraph(C), gamma = 1, n_runs = 5) cp$n_communities
Compares a subject's EMG coordination pattern against a reference (healthy baseline or normative data) to produce a 0-1 quality score.
mskCoordinationQualityScore( emg, reference_emg = NULL, hg = NULL, method = c("synergy_distance", "correlation_profile", "network_similarity") )mskCoordinationQualityScore( emg, reference_emg = NULL, hg = NULL, method = c("synergy_distance", "correlation_profile", "network_similarity") )
emg |
EMG data: matrix (time x channels), SummarizedExperiment, or vector. |
reference_emg |
Reference EMG matrix (or NULL for within-subject reference). |
hg |
An MSKHypergraph object (NULL loads default). |
method |
Character, comparison method: "synergy_distance" (default), "correlation_profile", or "network_similarity". |
A list with: quality_score (0-1, 1=perfect match to reference), component_scores, muscle_contributions.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_ref <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_ref) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") cqs <- mskCoordinationQualityScore(emg, emg_ref, hg = hg, method = "correlation_profile") cqsC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_ref <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_ref) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") cqs <- mskCoordinationQualityScore(emg, emg_ref, hg = hg, method = "correlation_profile") cqs
Detects compensatory movement patterns by comparing current EMG activation against a baseline, using MSK network topology to identify muscles that are structurally positioned to compensate for injured muscles.
mskDetectCompensation( emg, emg_baseline, injured_muscles, hg = NULL, emg_mapping = NULL, sr = NULL, z_threshold = 1.96, neighborhood_order = 2L )mskDetectCompensation( emg, emg_baseline, injured_muscles, hg = NULL, emg_mapping = NULL, sr = NULL, z_threshold = 1.96, neighborhood_order = 2L )
emg |
Current EMG data (matrix, SummarizedExperiment, or numeric vector). |
emg_baseline |
Baseline/pre-injury EMG data (same format as |
injured_muscles |
Character vector of injured muscle names or integer indices. |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
sr |
Optional sampling rate override. |
z_threshold |
Numeric, z-score threshold for flagging compensation (default: 1.96). |
neighborhood_order |
Integer, MSK graph distance to search for compensators (default: 2). |
An S3 object of class "MSKCompensation" with:
Data.frame of muscles showing compensatory activation
Data.frame of injured muscle activation status
Data.frame of neighbors that did NOT compensate
Proportion of neighbors showing compensation
List with neighborhood info
Map EMG channels to MSK muscles via emgToMSKMapping
Compute RMS activation for both current and baseline
Compute z-score of change: z = (current_rms - baseline_rms) / baseline_sd
Identify MSK neighbors of injured muscles within neighborhood_order hops
Flag muscles where they are neighbors AND z-score > z_threshold
Also detect decreased activation in injured muscles (z < -z_threshold)
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- mskDetectCompensation(emg, emg_base, "Biceps Brachii", hg = hg) print(result)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) set.seed(1); emg_base <- matrix(abs(rnorm(300 * 4)), 300, 4) set.seed(2); emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- colnames(emg_base) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- mskDetectCompensation(emg, emg_base, "Biceps Brachii", hg = hg) print(result)
Identifies when recovery has plateaued using rolling window slope analysis or change rate assessment.
mskDetectRecoveryPlateau( tracker, window = 3L, min_slope = NULL, method = c("slope", "change_rate") )mskDetectRecoveryPlateau( tracker, window = 3L, min_slope = NULL, method = c("slope", "change_rate") )
tracker |
An |
window |
Integer, number of consecutive timepoints to assess (default: 3). |
min_slope |
Numeric, minimum slope to be considered "improving" (default: NULL, auto-computed from data variability). |
method |
Character, "slope" (default) or "change_rate". |
An S3 object of class "MSKRecoveryPlateau" with:
logical per muscle
timepoint index where plateau begins (NA if none)
data.frame with per-window slopes
character ("continue", "modify_protocol", "reassess")
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) plateau <- mskDetectRecoveryPlateau(tracker, window = 2) plateauC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) plateau <- mskDetectRecoveryPlateau(tracker, window = 2) plateau
Compares first and last timepoint to determine if each muscle shows clinically meaningful change beyond the MDC threshold.
mskDetectResponderStatus( tracker, mdc_result = NULL, threshold_type = c("mdc", "effect_size"), instrument = NULL, population = NULL, direction = c("increase", "decrease") )mskDetectResponderStatus( tracker, mdc_result = NULL, threshold_type = c("mdc", "effect_size"), instrument = NULL, population = NULL, direction = c("increase", "decrease") )
tracker |
An |
mdc_result |
An |
threshold_type |
Character, "mdc" (default) or "effect_size" (Cohen's d > 0.8). |
instrument |
Optional validated instrument id (e.g. |
population |
Optional population stratum for the instrument's clinimetric
lookup (used only with |
direction |
|
When a validated instrument and population are supplied and the
PhysioClinical package is available, classification is delegated to
PhysioClinical::classifyResponder(), which applies the published dual
MDC-vs-MCID rule (adding a clinical_class column). Otherwise the
single-threshold path (MDC or effect size) is used.
An S3 object of class "MSKResponderStatus" with:
data.frame (muscle, status, change, threshold, effect_size_d; plus clinical_class when delegated)
counts of responders/non-responders/deteriorated
"responder" if majority of muscles improve
"clinimetric" (delegated) or "threshold"
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) status <- mskDetectResponderStatus(tracker) statusC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) status <- mskDetectResponderStatus(tracker) status
Tests whether a set of muscles is enriched for specific anatomical or functional annotations (actions, nerves, body regions, spinal levels) compared to the full background set of muscles in the hypergraph.
mskEnrichKG( muscles, annotation_type = c("action", "nerve", "body_region", "spinal_level"), hub = NULL, hg = NULL )mskEnrichKG( muscles, annotation_type = c("action", "nerve", "body_region", "spinal_level"), hub = NULL, hg = NULL )
muscles |
Character or integer vector identifying muscles to test. |
annotation_type |
Character, one of |
hub |
A PhysioAnnotationHub object (NULL loads default). |
hg |
An MSKHypergraph object (NULL loads default). |
A data.frame with columns:
Annotation term
Number of query muscles with this term
Number of background muscles with this term
Expected count under null
Ratio of observed to expected
P-value from hypergeometric test
Logical, TRUE if p < 0.05
Enrichment analysis depends on the completeness and accuracy of the underlying knowledge graph annotations. P-values are not adjusted for multiple testing; consider applying Bonferroni or FDR correction when testing multiple annotation types simultaneously.
if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { enrichment <- mskEnrichKG(c("Biceps Brachii", "Deltoid", "Trapezius"), annotation_type = "nerve") enrichment[enrichment$significant, ] }if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { enrichment <- mskEnrichKG(c("Biceps Brachii", "Deltoid", "Trapezius"), annotation_type = "nerve") enrichment[enrichment$significant, ] }
Takes the output of mskPredictFunctionalOutcome and generates
time-based milestones for tracking rehabilitation progress.
mskFunctionalMilestones( outcome_prediction, n_milestones = 4L, milestone_type = c("clinical", "linear", "accelerating") )mskFunctionalMilestones( outcome_prediction, n_milestones = 4L, milestone_type = c("clinical", "linear", "accelerating") )
outcome_prediction |
An |
n_milestones |
Integer, number of milestones to generate (default: 4). |
milestone_type |
Character: "linear" (evenly spaced), "accelerating" (front-loaded), or "clinical" (based on standard rehab phases). |
An S3 object of class "MSKFunctionalMilestones" with:
data.frame with milestone_id, week, target_rom, target_strength, target_function, phase_name, description
total timeline in weeks
character
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) milestones <- mskFunctionalMilestones(outcome, milestone_type = "clinical") print(milestones)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) milestones <- mskFunctionalMilestones(outcome, milestone_type = "clinical") print(milestones)
Tests the correspondence between MSK network community structure and the motor cortex homunculus. Reproduces Fig 4b results. Target: F(1,19) = 21.3, R² = 0.52, p < 0.001
mskHomuncCorrelation(membership = NULL, homunculus_data = NULL)mskHomuncCorrelation(membership = NULL, homunculus_data = NULL)
membership |
Named integer vector of community assignments. |
homunculus_data |
Optional data.frame with columns: homunc_area, dev_ratio. If NULL, loads the built-in validation data. |
A list with regression result and deviation ratio by category.
Murphy AC et al. (2018) PLOS Biology.
res <- mskHomuncCorrelation() res$r_squaredres <- mskHomuncCorrelation() res$r_squared
Constructs a musculoskeletal hypergraph from an incidence matrix. In this hypergraph, bones are vertices and muscles are hyperedges. A muscle (hyperedge) connects all bones it attaches to.
MSKHypergraph(C = NULL, muscle_meta = NULL)MSKHypergraph(C = NULL, muscle_meta = NULL)
C |
A matrix or sparse Matrix (bones x muscles) where C[i,j]=1 means muscle j attaches to bone i. If NULL, loads the built-in data. |
muscle_meta |
Optional data.frame with muscle metadata. |
An S3 object of class "MSKHypergraph" containing:
Sparse incidence matrix (bones x muscles)
Number of bones (vertices)
Number of muscles (hyperedges)
Character vector of bone names
Character vector of muscle names
Muscle metadata if provided
Murphy AC et al. (2018) PLOS Biology 16(1): e2002811.
## Tiny synthetic incidence matrix: 4 bones x 4 muscles (a ring) C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) hg ## Or load the bundled Murphy et al. (2018) dataset hg_full <- MSKHypergraph() hg_full$n_bones## Tiny synthetic incidence matrix: 4 bones x 4 muscles (a ring) C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) hg ## Or load the bundled Murphy et al. (2018) dataset hg_full <- MSKHypergraph() hg_full$n_bones
Computes the impact deviation for each muscle, which is the difference between the observed impact score and the expected impact score for a muscle of that degree, expressed in standard deviations.
mskImpactDeviation(impact_scores, hg, null_scores = NULL)mskImpactDeviation(impact_scores, hg, null_scores = NULL)
impact_scores |
Named numeric vector of impact scores (from mskImpactScoreAll). |
hg |
An MSKHypergraph object. |
null_scores |
Optional matrix of null model impact scores (muscles x null_runs). If NULL, deviation is computed relative to a degree-based regression. |
The expected impact score is estimated from null model simulations or from a regression of impact score vs. degree.
A named numeric vector of impact deviations.
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) scores <- mskImpactScoreAll(mskSimulate(hg, n_steps = 50), verbose = FALSE) mskImpactDeviation(scores, hg)C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) scores <- mskImpactScoreAll(mskSimulate(hg, n_steps = 50), verbose = FALSE) mskImpactDeviation(scores, hg)
Reproduces the key result from Fig 3b: correlation between impact deviation and clinical muscle injury recovery time. Target: F(1,12) = 37.3, R² = 0.757, p < 0.0001
mskImpactRecoveryModel(impact_deviation = NULL, recovery_data = NULL)mskImpactRecoveryModel(impact_deviation = NULL, recovery_data = NULL)
impact_deviation |
Named numeric vector of impact deviations. |
recovery_data |
Optional data.frame with columns: recovery_time, impact_deviation, weight. If NULL, loads the built-in validation data. |
Result from mskRobustRegression with additional paper comparison.
Murphy AC et al. (2018) PLOS Biology Table 4.
res <- mskImpactRecoveryModel() res$r_squaredres <- mskImpactRecoveryModel() res$r_squared
Perturbs a single muscle and simulates the network dynamics to compute the total displacement of all bones (impact score). The perturbation is applied in the 4th spatial dimension to avoid directional artifacts.
mskImpactScore(sim, muscle_index)mskImpactScore(sim, muscle_index)
sim |
An MSKSimulation object. |
muscle_index |
Integer index of the muscle to perturb. |
Numeric impact score (total displacement summed over all bones).
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) sim <- mskSimulate(MSKHypergraph(C), n_steps = 50) mskImpactScore(sim, muscle_index = 1)C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) sim <- mskSimulate(MSKHypergraph(C), n_steps = 50) mskImpactScore(sim, muscle_index = 1)
Runs perturbation analysis for all 270 muscles and returns impact scores. This is the main computational function for reproducing Fig. 3a of the paper.
mskImpactScoreAll(sim = NULL, verbose = TRUE)mskImpactScoreAll(sim = NULL, verbose = TRUE)
sim |
An MSKSimulation object. If NULL, creates one with default parameters. |
verbose |
Logical, print progress (default: TRUE). |
A named numeric vector of impact scores for each muscle.
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) sim <- mskSimulate(MSKHypergraph(C), n_steps = 50) mskImpactScoreAll(sim, verbose = FALSE)C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) sim <- mskSimulate(MSKHypergraph(C), n_steps = 50) mskImpactScoreAll(sim, verbose = FALSE)
Computes a personalized injury risk profile for all muscles based on MSK network topology and patient characteristics.
mskInjuryRiskProfile(patient_data, hg = NULL, sim = NULL)mskInjuryRiskProfile(patient_data, hg = NULL, sim = NULL)
patient_data |
A list with patient characteristics:
|
hg |
An MSKHypergraph object (NULL loads default). |
sim |
An MSKSimulation object (NULL creates default). |
An S3 object of class "MSKInjuryRiskProfile" with:
Data frame with muscle name, base risk, adjusted risk
Input patient data
Applied adjustment factors
Patient factor adjustments (age, BMI, activity level) are heuristic multipliers, not derived from validated epidemiological models. This is a research exploration tool, not a clinical diagnostic.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) profile <- mskInjuryRiskProfile(list(age = 45, activity_level = "active"), hg = hg) profileC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) profile <- mskInjuryRiskProfile(list(age = 45, activity_level = "active"), hg = hg) profile
Produces a comprehensive summary of the MSK network annotated with knowledge graph data. Profiles all communities, identifies cross-community nerve pathways, and summarizes annotation coverage.
mskKGSummary(hg = NULL, hub = NULL, gamma = 4.3)mskKGSummary(hg = NULL, hub = NULL, gamma = 4.3)
hg |
An MSKHypergraph object (NULL loads default). |
hub |
A PhysioAnnotationHub object (NULL loads default). |
gamma |
Numeric, resolution parameter for community detection (default: 4.3). |
An S3 object of class "MSKKGSummary" with:
Annotated MSKHypergraph
List of MSKCommunityProfile objects
Number of detected communities
Data frame of nerves spanning communities
Annotation match statistics
This is a summary tool combining network topology with knowledge graph annotations. Community boundaries and annotation mappings are model-derived. Use as a research exploration tool, not for clinical decision-making.
if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { s <- mskKGSummary() print(s) }if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { s <- mskKGSummary() print(s) }
Maps an MSK-network functional-outcome prediction to a validated clinical
outcome measure (COM) instrument and its WHO ICF category tags (via
PhysioAnnotationHub::tagICF(), when that package is installed). This
bridges an abstract network prediction to a documented, ICF-anchored clinical
instrument.
mskLinkOutcomeMeasure( outcome, instrument_id, measure = c("function", "rom", "strength"), hub = NULL )mskLinkOutcomeMeasure( outcome, instrument_id, measure = c("function", "rom", "strength"), hub = NULL )
outcome |
An |
instrument_id |
A clinical outcome-measure instrument id (e.g.
|
measure |
Which aggregate prediction to link when |
hub |
Optional |
A data.frame with instrument_id, measure,
predicted_value, and icf_code (one row per ICF tag; a single
NA icf_code when no ICF link is available).
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) mskLinkOutcomeMeasure(outcome, "fma_ue", measure = "strength")C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) mskLinkOutcomeMeasure(outcome, "fma_ue", measure = "strength")
Computes per-timepoint RMS activation, synergy decomposition (W/H/VAF), and optionally CMC for a series of EMG measurements taken at different timepoints during rehabilitation.
mskLongitudinalTracker( timepoints, hg = NULL, emg_mapping = NULL, metrics = c("rms", "synergy"), sr = NULL )mskLongitudinalTracker( timepoints, hg = NULL, emg_mapping = NULL, metrics = c("rms", "synergy"), sr = NULL )
timepoints |
Named list of EMG matrices (time x channels). Names should be timepoint labels (e.g., "T0", "T1", "T2"). |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
metrics |
Character vector of metrics to compute (default: c("rms", "synergy")). Options: "rms", "synergy", "mean_activation", "peak_activation". |
sr |
Optional sampling rate in Hz (overrides detected values). |
An S3 object of class "MSKLongitudinalTracker" with:
data.frame (timepoint, muscle, metric_name, value)
character vector
integer
list of synergy results per timepoint
matrix (n_muscles x n_timepoints)
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) print(tracker)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) print(tracker)
Computes ICC (intraclass correlation coefficient) across timepoints for each metric, then derives SEM and MDC values for clinimetric assessment.
mskMinimalDetectableChange( tracker, method = "standard", confidence = 0.95, icc_method = "ICC(2,1)" )mskMinimalDetectableChange( tracker, method = "standard", confidence = 0.95, icc_method = "ICC(2,1)" )
tracker |
An |
method |
Character, method for MDC computation (default: "standard"). |
confidence |
Numeric, confidence level for MDC (default: 0.95). |
icc_method |
Character, ICC formula variant (default: "ICC(2,1)"). |
An S3 object of class "MSKMinimalDetectableChange" with:
data.frame (muscle, metric, ICC, SEM, MDC, pooled_SD)
numeric
character
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) mdc <- mskMinimalDetectableChange(tracker, confidence = 0.95) mdcC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) mdc <- mskMinimalDetectableChange(tracker, confidence = 0.95) mdc
Computes the modularity Q of a partition with resolution parameter. Q = (1/2m) * sum_ij (A_ij - gamma * k_i * k_j / (2m)) * delta(c_i, c_j)
mskModularity(A, membership, gamma = 1)mskModularity(A, membership, gamma = 1)
A |
Adjacency matrix. |
membership |
Integer vector of community assignments. |
gamma |
Resolution parameter (default: 1.0). |
Numeric modularity value.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) A <- projectMuscleGraph(hg) cm <- mskCommunityDetect(hg, gamma = 1) mskModularity(A, cm$membership, gamma = 1)C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) A <- projectMuscleGraph(hg) cm <- mskCommunityDetect(hg, gamma = 1) mskModularity(A, cm$membership, gamma = 1)
Computes standard graph metrics for the projected bone or muscle graph.
mskNetworkMetrics(hg, type = c("bone", "muscle"))mskNetworkMetrics(hg, type = c("bone", "muscle"))
hg |
An MSKHypergraph object. |
type |
Character, "bone" or "muscle" projection. |
A list with: degree, strength, clustering_coef, density
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) m <- mskNetworkMetrics(MSKHypergraph(C), type = "bone") m$density m$degreeC <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) m <- mskNetworkMetrics(MSKHypergraph(C), type = "bone") m$density m$degree
Compares corticomuscular coherence (CMC) and directional coupling between two timepoints to quantify neural adaptation during rehabilitation.
mskNeuralAdaptationIndex( cmc_t0, cmc_t1, directional_t0 = NULL, directional_t1 = NULL )mskNeuralAdaptationIndex( cmc_t0, cmc_t1, directional_t0 = NULL, directional_t1 = NULL )
cmc_t0 |
CMC result at baseline (MSKNeuromechCMC object or coherence matrix). |
cmc_t1 |
CMC result at follow-up. |
directional_t0 |
Directional coupling at baseline (MSKNeuromechDirectional object or NULL). |
directional_t1 |
Directional coupling at follow-up (or NULL). |
A list with: cmc_change, directional_change, adaptation_index, interpretation ("improving"/"stable"/"declining").
set.seed(1) cmc_t0 <- matrix(abs(rnorm(9)), 3, 3) cmc_t1 <- matrix(abs(rnorm(9)), 3, 3) nai <- mskNeuralAdaptationIndex(cmc_t0, cmc_t1) naiset.seed(1) cmc_t0 <- matrix(abs(rnorm(9)), 3, 3) cmc_t1 <- matrix(abs(rnorm(9)), 3, 3) nai <- mskNeuralAdaptationIndex(cmc_t0, cmc_t1) nai
Creates multiple null hypergraphs and computes impact scores for each. This is used to compute impact deviations relative to the null distribution.
mskNullEnsemble(hg = NULL, n_null = 100L, sim_params = list(), verbose = TRUE)mskNullEnsemble(hg = NULL, n_null = 100L, sim_params = list(), verbose = TRUE)
hg |
An MSKHypergraph object. |
n_null |
Number of null models to generate (default: 100). |
sim_params |
List of simulation parameters (dt, n_steps, beta). |
verbose |
Logical, print progress. |
A list with:
Matrix of impact scores (n_muscles x n_null)
List, impact scores grouped by degree for each null
Number of null models
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) set.seed(1) ens <- mskNullEnsemble(MSKHypergraph(C), n_null = 3, sim_params = list(n_steps = 50), verbose = FALSE) dim(ens$null_scores)C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) set.seed(1) ens <- mskNullEnsemble(MSKHypergraph(C), n_null = 3, sim_params = list(n_steps = 50), verbose = FALSE) dim(ens$null_scores)
Creates a randomized version of the musculoskeletal hypergraph by rewiring muscle-bone connections while preserving each muscle's degree (number of bones it connects to). This is done by randomly reassigning which bones each muscle attaches to, within anatomical categories if specified.
mskNullHypergraph(hg = NULL, preserve_category = FALSE)mskNullHypergraph(hg = NULL, preserve_category = FALSE)
hg |
An MSKHypergraph object. If NULL, loads built-in data. |
preserve_category |
Logical, whether to preserve anatomical category during rewiring (default: FALSE for the basic model). |
An MSKHypergraph object with rewired connections.
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) set.seed(1) null_hg <- mskNullHypergraph(MSKHypergraph(C)) null_hg$n_musclesC <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) set.seed(1) null_hg <- mskNullHypergraph(MSKHypergraph(C)) null_hg$n_muscles
Computes confidence intervals for functional outcome predictions using bootstrap resampling or analytical delta method.
mskOutcomeConfidenceInterval( predictions, method = c("bootstrap", "analytical"), n_boot = 100L, confidence_level = 0.95 )mskOutcomeConfidenceInterval( predictions, method = c("bootstrap", "analytical"), n_boot = 100L, confidence_level = 0.95 )
predictions |
An |
method |
Character: "bootstrap" (default) or "analytical". |
n_boot |
Integer, number of bootstrap resamples (default: 100). |
confidence_level |
Numeric, confidence level (default: 0.95). |
A list with:
data.frame with outcome, point_estimate, lower, upper, width
character
integer (for bootstrap)
numeric, mean CI width
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) ci <- mskOutcomeConfidenceInterval(outcome) ciC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) ci <- mskOutcomeConfidenceInterval(outcome) ci
Orchestrator function that combines functional outcome prediction, milestones, and optional reassessment into a comprehensive report.
mskOutcomeReport(outcome, milestones = NULL, reassessment = NULL, hg = NULL)mskOutcomeReport(outcome, milestones = NULL, reassessment = NULL, hg = NULL)
outcome |
An |
milestones |
Optional |
reassessment |
Optional |
hg |
An MSKHypergraph object (NULL loads default). |
An S3 object of class "MSKOutcomeReport" with all sub-results.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) milestones <- mskFunctionalMilestones(outcome) report <- mskOutcomeReport(outcome, milestones) print(report)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) milestones <- mskFunctionalMilestones(outcome) report <- mskOutcomeReport(outcome, milestones) print(report)
Orchestrator function that calls all clinical bridge functions with shared hypergraph and simulation objects, producing a comprehensive report.
mskOutcomeSummary(injury_muscles, patient_data = NULL, hg = NULL, sim = NULL)mskOutcomeSummary(injury_muscles, patient_data = NULL, hg = NULL, sim = NULL)
injury_muscles |
Character or integer vector identifying injured muscles. |
patient_data |
Optional list with patient characteristics (see
|
hg |
An MSKHypergraph object (NULL loads default). |
sim |
An MSKSimulation object (NULL creates default). |
An S3 object of class "MSKOutcomeSummary" containing:
MSKClinicalPrediction object
Recovery timeline data.frame
MSKInjuryRiskProfile object (if patient_data provided)
MSKRehabProtocol object
All predictions are model-based estimates using MSK network topology. Recovery model was validated on 14 aggregate muscle groups. Patient factors are heuristic. This is a research exploration tool, not a clinical diagnostic.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) s <- mskOutcomeSummary("Biceps Brachii", hg = hg) print(s)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) s <- mskOutcomeSummary("Biceps Brachii", hg = hg) print(s)
Finds the shortest anatomical pathway between two entities (muscles, bones, nerves, etc.) in the knowledge graph. Useful for tracing innervation chains, biomechanical linkages, and anatomical relationships.
mskPathwayQuery(from, to, hub = NULL, max_depth = 5L)mskPathwayQuery(from, to, hub = NULL, max_depth = 5L)
from |
Character, the source entity name (e.g., "Biceps Brachii"). |
to |
Character, the target entity name (e.g., "C5"). |
hub |
A PhysioAnnotationHub object (NULL loads default). |
max_depth |
Integer, maximum BFS depth (default: 5). |
A list with:
Character vector of entities along the path
Character vector of relationship types between entities
Integer, path length
Human-readable path description
Logical, whether a path was found
Pathways reflect relationships encoded in the knowledge graph. The shortest path in the KG may not correspond to the most clinically relevant connection. Always verify pathway interpretations against anatomical references.
if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { path <- mskPathwayQuery("Biceps Brachii", "C5") cat(path$description, "\n") }if (requireNamespace("PhysioAnnotationHub", quietly = TRUE)) { path <- mskPathwayQuery("Biceps Brachii", "C5") cat(path$description, "\n") }
Predicts clinical functional outcomes (ROM, strength, composite function score) for injured muscles based on MSK network topology, impact analysis, and optional EMG activation data. Provides evidence-based regression models with confidence intervals.
mskPredictFunctionalOutcome( injured_muscles, hg = NULL, outcome_type = c("all", "rom", "strength", "function"), patient_factors = NULL, emg = NULL, emg_mapping = NULL, confidence_level = 0.95 )mskPredictFunctionalOutcome( injured_muscles, hg = NULL, outcome_type = c("all", "rom", "strength", "function"), patient_factors = NULL, emg = NULL, emg_mapping = NULL, confidence_level = 0.95 )
injured_muscles |
Character or integer vector identifying injured muscles. |
hg |
An MSKHypergraph object (NULL loads default 173-bone/270-muscle network). |
outcome_type |
Character: "rom" (range of motion), "strength", "function" (composite functional score), or "all" (default). |
patient_factors |
Optional list with: age (numeric), sex (character), bmi (numeric), activity_level (character), injury_severity ("mild"/"moderate"/"severe"). |
emg |
Optional EMG data (SummarizedExperiment, matrix, or numeric vector) for activation-based prediction refinement. |
emg_mapping |
Optional pre-computed data.frame from |
confidence_level |
Numeric, confidence level for CIs (default: 0.95). |
An S3 object of class "MSKFunctionalOutcome" with:
data.frame with muscle, outcome_type, predicted_value, lower_ci, upper_ci, unit
list with overall_rom, overall_strength, overall_function
estimated weeks to reach 90 percent of predicted outcome
numeric
character
list
The prediction models are based on MSK network topology and heuristic adjustments. They are not independently validated clinical tools. Use as a research exploration tool, not a clinical diagnostic.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) print(outcome)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) print(outcome)
Compares current physiological measurements to a previous assessment, computing progress metrics and generating recommendations.
mskReassess( current_data, previous_assessment, hg = NULL, emg_mapping = NULL, sr = NULL )mskReassess( current_data, previous_assessment, hg = NULL, emg_mapping = NULL, sr = NULL )
current_data |
A list with: emg (required, matrix or SummarizedExperiment), and optional eeg, kinematics, force. |
previous_assessment |
A previous |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed EMG-to-MSK mapping. |
sr |
Optional numeric sampling rate. |
An S3 object of class "MSKReassessment" with:
data.frame of current activation/synergy metrics
data.frame from previous assessment
data.frame with metric, previous, current, change, pct_change, significant
character vector per metric
weighted average progress score (0-100)
character
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") reassessment <- mskReassess(list(emg = emg), outcome, hg = hg) print(reassessment)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) outcome <- mskPredictFunctionalOutcome("Biceps Brachii", hg = hg) set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") reassessment <- mskReassess(list(emg = emg), outcome, hg = hg) print(reassessment)
Generates a week-by-week recovery timeline based on exponential decay of network impact, with clinical phase assignments.
mskRecoveryTimeline(injury_muscles, hg = NULL, sim = NULL, n_weeks = 12L)mskRecoveryTimeline(injury_muscles, hg = NULL, sim = NULL, n_weeks = 12L)
injury_muscles |
Character or integer vector identifying injured muscles. |
hg |
An MSKHypergraph object (NULL loads default). |
sim |
An MSKSimulation object (NULL creates default). |
n_weeks |
Integer, number of weeks to project (default: 12). |
A data.frame with columns: week, muscle, remaining_impact_pct, phase, milestone.
Phase assignments are based on exponential decay of network impact scores, not empirical clinical data. They should be interpreted as model-based estimates for research purposes only.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) timeline <- mskRecoveryTimeline("Biceps Brachii", hg = hg) timelineC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) timeline <- mskRecoveryTimeline("Biceps Brachii", hg = hg) timeline
Fits parametric recovery curves (exponential, sigmoid, or linear) to longitudinal muscle activation data from an MSKLongitudinalTracker.
mskRecoveryTrajectoryFit( tracker, model = c("exponential", "sigmoid", "linear"), muscle_subset = NULL, partial_pool = FALSE )mskRecoveryTrajectoryFit( tracker, model = c("exponential", "sigmoid", "linear"), muscle_subset = NULL, partial_pool = FALSE )
tracker |
An |
model |
Character, recovery curve model: "exponential" (default), "sigmoid", or "linear". |
muscle_subset |
Optional character vector of muscle names to fit (NULL = all muscles). |
partial_pool |
Logical; when |
An S3 object of class "MSKRecoveryTrajectory" with:
list per muscle with coefficients, residuals, R-squared, AIC
matrix of predicted values (muscles x timepoints)
character
numeric vector (slope at midpoint for each muscle)
estimated time to 90 percent of asymptotic recovery
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) traj <- mskRecoveryTrajectoryFit(tracker, model = "linear") trajC <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) mk <- function(s) { set.seed(s); m <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(m) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii"); m } tracker <- mskLongitudinalTracker(list(T0 = mk(1), T1 = mk(2), T2 = mk(3)), hg = hg) traj <- mskRecoveryTrajectoryFit(tracker, model = "linear") traj
Generates a phased rehabilitation protocol using MSK network community structure and shortest path distances to determine exercise progression.
mskRehabProtocol(injury_muscles, hg = NULL, n_phases = 3L)mskRehabProtocol(injury_muscles, hg = NULL, n_phases = 3L)
injury_muscles |
Character or integer vector identifying injured muscles. |
hg |
An MSKHypergraph object (NULL loads default). |
n_phases |
Integer, number of rehabilitation phases (default: 3). |
An S3 object of class "MSKRehabProtocol" with per-phase muscle lists.
Phase ordering is based on network topology (community membership and shortest path distance), not validated rehabilitation protocols. Use as a research exploration tool, not clinical guidance.
C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) protocol <- mskRehabProtocol("Biceps Brachii", hg = hg) print(protocol)C <- matrix(c(1,1,0,0,0, 0,1,1,0,0, 0,0,1,1,0, 0,0,0,1,1, 1,0,0,0,1, 0,1,0,1,0), nrow = 6, byrow = TRUE, dimnames = list(paste0("bone", 1:6), c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii", "Brachialis"))) hg <- MSKHypergraph(C) protocol <- mskRehabProtocol("Biceps Brachii", hg = hg) print(protocol)
Master function that runs the complete analysis pipeline and compares results to the paper's reported values. This is the main validation function.
mskReproducePaper(run_simulation = FALSE, verbose = TRUE)mskReproducePaper(run_simulation = FALSE, verbose = TRUE)
run_simulation |
Logical, run the full simulation (slow, default: FALSE). If FALSE, uses validation data to reproduce statistical analyses. |
verbose |
Logical, print progress and comparison. |
A list with all results and paper comparisons.
res <- mskReproducePaper(verbose = FALSE) res$recovery$r_squaredres <- mskReproducePaper(verbose = FALSE) res$recovery$r_squared
Fits a robust linear model (using MASS::rlm if available, otherwise lm). Supports weighted regression as used in the paper.
mskRobustRegression(x, y, weights = NULL)mskRobustRegression(x, y, weights = NULL)
x |
Numeric vector of predictor values. |
y |
Numeric vector of response values. |
weights |
Optional numeric vector of weights. |
A list with:
Named vector (intercept, slope)
R-squared value
F-statistic
p-value for the regression
Residual values
Fitted values
The fitted model object
set.seed(1) x <- 1:20 y <- 2 * x + rnorm(20, sd = 2) fit <- mskRobustRegression(x, y) fit$coefficients fit$r_squaredset.seed(1) x <- 1:20 y <- 2 * x + rnorm(20, sd = 2) fit <- mskRobustRegression(x, y) fit$coefficients fit$r_squared
Computes shortest path distance matrix for the projected graph.
mskShortestPaths(hg, type = c("bone", "muscle"))mskShortestPaths(hg, type = c("bone", "muscle"))
hg |
An MSKHypergraph object. |
type |
Character, "bone" or "muscle" projection. |
A numeric matrix of shortest path distances.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) mskShortestPaths(MSKHypergraph(C), type = "bone")C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) mskShortestPaths(MSKHypergraph(C), type = "bone")
Runs a damped harmonic oscillator simulation on the musculoskeletal network. Each muscle is modeled as a spring connecting its attached bones, and each bone is a unit-mass point particle. The system is evolved under perturbation and the total displacement is computed as the impact score.
mskSimulate( hg = NULL, dt = 0.01, n_steps = 500L, beta = 1, perturbation_magnitude = 1 )mskSimulate( hg = NULL, dt = 0.01, n_steps = 500L, beta = 1, perturbation_magnitude = 1 )
hg |
An MSKHypergraph object. If NULL, loads the built-in data. |
dt |
Time step for integration (default: 0.01). |
n_steps |
Number of time steps to integrate (default: 500). |
beta |
Damping coefficient (default: 1.0). |
perturbation_magnitude |
Magnitude of the 4th-dimension perturbation (default: 1.0). |
An S3 object of class "MSKSimulation" with:
The MSKHypergraph used
Named vector of spring constants per muscle
List of simulation parameters
Murphy AC et al. (2018) PLOS Biology 16(1): e2002811.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) sim <- mskSimulate(MSKHypergraph(C), n_steps = 50) simC <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) sim <- mskSimulate(MSKHypergraph(C), n_steps = 50) sim
Compares two synergy decompositions by aligning synergy weight vectors and computing a change index reflecting structural reorganization.
mskSynergyChangeIndex( synergy_t0, synergy_t1, method = c("cosine", "correlation", "procrustes") )mskSynergyChangeIndex( synergy_t0, synergy_t1, method = c("cosine", "correlation", "procrustes") )
synergy_t0 |
An |
synergy_t1 |
An |
method |
Character, comparison method: "cosine" (default), "correlation", or "procrustes". |
A list with: global_change_index (0=identical, 1=maximally different), per_synergy_change, alignment (permutation used).
set.seed(1) synergy_t0 <- matrix(abs(rnorm(20)), 5, 4) synergy_t1 <- matrix(abs(rnorm(20)), 5, 4) sci <- mskSynergyChangeIndex(synergy_t0, synergy_t1, method = "cosine") sciset.seed(1) synergy_t0 <- matrix(abs(rnorm(20)), 5, 4) synergy_t1 <- matrix(abs(rnorm(20)), 5, 4) sci <- mskSynergyChangeIndex(synergy_t0, synergy_t1, method = "cosine") sci
Computes the z-scored Rand index between two partitions. Used to compare community structure with homunculus categories.
mskZRand(partition1, partition2)mskZRand(partition1, partition2)
partition1 |
Integer vector of community assignments. |
partition2 |
Integer vector of category assignments. |
Numeric z-Rand score. Values > 1.96 indicate significant similarity.
Traud et al. (2011) Physical Review E.
p1 <- c(1, 1, 1, 2, 2, 2) p2 <- c(1, 1, 2, 2, 2, 2) mskZRand(p1, p2)p1 <- c(1, 1, 1, 2, 2, 2) p2 <- c(1, 1, 2, 2, 2, 2) mskZRand(p1, p2)
Describes muscle-tendon-unit (MTU) length as a polynomial in the joint angles the muscle crosses:
Moment arms follow as -dL/dtheta_j (Menegaldo et al. 2004). Supply hill
parameters to enable normalized fiber length/velocity in
muscleTendonKinematics().
musclePath(joints, coefficients, slack_length = 0, hill = NULL, muscle = NULL)musclePath(joints, coefficients, slack_length = 0, hill = NULL, muscle = NULL)
joints |
Character vector of joint names the muscle crosses. |
coefficients |
Named list (by joint) of polynomial coefficient vectors
|
slack_length |
MTU length (m) at all-zero joint angles (default 0). |
hill |
Optional list with |
muscle |
Optional muscle name (metadata). |
A muscle_path object.
Menegaldo LL, et al. (2004). J Biomech 37(9):1447-1453.
muscleTendonKinematics(), defaultMusclePath()
# gastrocnemius crosses knee and ankle musclePath(c("knee", "ankle"), coefficients = list(knee = c(0.03), ankle = c(-0.05, 0.01)), slack_length = 0.42)# gastrocnemius crosses knee and ankle musclePath(c("knee", "ankle"), coefficients = list(knee = c(0.03), ankle = c(-0.05, 0.01)), slack_length = 0.42)
Evaluates a musclePath() over one or more poses: MTU length L(theta),
moment arms -dL/dtheta_j, and MTU lengthening velocity
dL/dt = -sum_j r_j * omega_j. When the model carries hill parameters,
rigid-tendon normalized fiber length and velocity are also returned.
muscleTendonKinematics( model, angles, angular_velocities = NULL, sampling_rate = NULL )muscleTendonKinematics( model, angles, angular_velocities = NULL, sampling_rate = NULL )
model |
A |
angles |
Joint angles in radians: a named numeric vector (single pose),
or an |
angular_velocities |
Joint angular velocities (rad/s), same shape as
|
sampling_rate |
Sampling rate (Hz) used to differentiate |
A muscle_tendon_kinematics object: a list with length (m),
velocity (m/s), moment_arms (n_frames x n_joints, m), joints, and,
when hill is set, norm_fiber_length and norm_fiber_velocity.
Zajac FE (1989). "Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control." Thelen DG (2003).
mp <- musclePath("ankle", coefficients = list(ankle = c(-0.05, 0.01)), slack_length = 0.3) muscleTendonKinematics(mp, angles = c(ankle = 0.1), angular_velocities = c(ankle = 2))mp <- musclePath("ankle", coefficients = list(ankle = c(-0.05, 0.01)), slack_length = 0.3) muscleTendonKinematics(mp, angles = c(ankle = 0.1), angular_velocities = c(ankle = 2))
Computes pairwise corticomuscular coherence (CMC) between EEG and EMG channels, maps EMG channels to MSK muscles, builds a CMC-based functional distance matrix, and compares it with the structural muscle adjacency via Mantel test.
neuromechCorticomuscularCoupling( eeg, emg, hg = NULL, freq_band = c(15, 35), eeg_channels = NULL, emg_mapping = NULL, sr_eeg = NULL, sr_emg = NULL, nperseg = 256L, n_perm = 999L )neuromechCorticomuscularCoupling( eeg, emg, hg = NULL, freq_band = c(15, 35), eeg_channels = NULL, emg_mapping = NULL, sr_eeg = NULL, sr_emg = NULL, nperseg = 256L, n_perm = 999L )
eeg |
EEG data: a SummarizedExperiment, numeric matrix (time x channels), or numeric vector. |
emg |
EMG data: a SummarizedExperiment, numeric matrix (time x channels), or numeric vector. |
hg |
An MSKHypergraph object (NULL loads default). |
freq_band |
Numeric vector of length 2, frequency band in Hz for CMC (default: c(15, 35) for beta range). |
eeg_channels |
Optional character vector of EEG channel names to use.
If NULL, motor cortex channels are auto-selected via |
emg_mapping |
Optional pre-computed data.frame from |
sr_eeg |
Optional sampling rate for EEG (overrides detected value). |
sr_emg |
Optional sampling rate for EMG (overrides detected value). |
nperseg |
Integer, segment length for Welch's method (default: 256). |
n_perm |
Integer, number of permutations for Mantel test (default: 999). |
An S3 object of class "MSKNeuromechCMC" with:
Coherence matrix (n_eeg x n_emg)
MSK muscle adjacency (matched subset)
Per-muscle mean CMC across EEG channels
Data.frame of significant EEG-EMG pairs
Mantel test result (correlation, p_value)
EMG-to-muscle mapping used
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") result <- neuromechCorticomuscularCoupling(eeg, emg, sr_eeg = 1000, sr_emg = 1000, nperseg = 128, n_perm = 19) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") result <- neuromechCorticomuscularCoupling(eeg, emg, sr_eeg = 1000, sr_emg = 1000, nperseg = 128, n_perm = 19) result
Computes directional (causal) connectivity between EEG and EMG channels using Granger causality or transfer entropy, distinguishing descending (cortical->muscle) from ascending (proprioceptive) pathways.
neuromechDirectionalCoupling( eeg, emg, hg = NULL, method = c("granger", "transfer_entropy", "both"), max_order_ms = 50, lag_ms = 20, n_bins = NULL, eeg_channels = NULL, emg_mapping = NULL, sr_eeg = NULL, sr_emg = NULL, n_perm = 999L, alpha = 0.05 )neuromechDirectionalCoupling( eeg, emg, hg = NULL, method = c("granger", "transfer_entropy", "both"), max_order_ms = 50, lag_ms = 20, n_bins = NULL, eeg_channels = NULL, emg_mapping = NULL, sr_eeg = NULL, sr_emg = NULL, n_perm = 999L, alpha = 0.05 )
eeg |
EEG data: SummarizedExperiment, matrix (time x channels), or vector. |
emg |
EMG data: SummarizedExperiment, matrix (time x channels), or vector. |
hg |
An MSKHypergraph object (NULL loads default). |
method |
Character: "granger" (default), "transfer_entropy", or "both". |
max_order_ms |
Numeric, maximum lag in ms for Granger (default: 50). |
lag_ms |
Numeric, TE lag in ms (default: 20). |
n_bins |
Integer, bins for TE discretization (NULL for auto). |
eeg_channels |
Optional character vector of EEG channels to use. |
emg_mapping |
Optional pre-computed data.frame from |
sr_eeg |
Optional sampling rate for EEG. |
sr_emg |
Optional sampling rate for EMG. |
n_perm |
Integer, permutations for Mantel test (default: 999). |
alpha |
Numeric, significance threshold (default: 0.05). |
An S3 object of class "MSKNeuromechDirectional" with:
Matrix (n_eeg x n_emg): EEG->EMG values
Matrix (n_emg x n_eeg): EMG->EEG values
descending - t(ascending)
Data.frame of significant descending pairs
Data.frame of significant ascending pairs
Per-muscle mean(descending)/mean(ascending)
Data.frame classifying each pair
Mantel test result
Method used
List of parameters used
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") result <- neuromechDirectionalCoupling(eeg, emg, method = "granger", sr_eeg = 1000, sr_emg = 1000, n_perm = 19) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") result <- neuromechDirectionalCoupling(eeg, emg, method = "granger", sr_eeg = 1000, sr_emg = 1000, n_perm = 19) result
Computes the electromechanical delay (EMD) between EMG and kinematic signals for muscle-bone pairs connected in the MSK hypergraph. Optionally correlates EMD with MSK network distance.
neuromechElectromechanicalDelay( emg, kinematics, hg = NULL, emg_mapping = NULL, kin_mapping = NULL, sr = NULL, max_lag_ms = 200, window_sec = NULL, n_perm = 999L )neuromechElectromechanicalDelay( emg, kinematics, hg = NULL, emg_mapping = NULL, kin_mapping = NULL, sr = NULL, max_lag_ms = 200, window_sec = NULL, n_perm = 999L )
emg |
EMG data: SummarizedExperiment, matrix (time x channels), or vector. |
kinematics |
Kinematic data: matrix (time x segments) or SummarizedExperiment. |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
kin_mapping |
Optional pre-computed data.frame from
|
sr |
Optional sampling rate (overrides detected value). |
max_lag_ms |
Numeric, maximum lag in milliseconds (default: 200). |
window_sec |
Optional numeric, window size in seconds for sliding window EMD analysis (NULL for global only). |
n_perm |
Integer, number of permutations for correlation test (default: 999). |
A list with:
Data.frame with muscle, bone, emd_ms, peak_correlation per pair
Corresponding MSK shortest path distance
Pearson correlation between EMD and network distance
Permutation p-value for the correlation
Data.frame with per-window EMD (if window_sec set)
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(4) kin <- matrix(cumsum(rnorm(300)) + sin(seq(0, 6, length.out = 300)), 300, 1) colnames(kin) <- "elbow" result <- neuromechElectromechanicalDelay(emg, kin, sr = 1000, n_perm = 19) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(4) kin <- matrix(cumsum(rnorm(300)) + sin(seq(0, 6, length.out = 300)), 300, 1) colnames(kin) <- "elbow" result <- neuromechElectromechanicalDelay(emg, kin, sr = 1000, n_perm = 19) result
Combines EMG activation, kinematic stress, and force data to compute a multi-source vulnerability score for each muscle. Propagates kinematic stress through the MSK incidence matrix and weights by impact deviation.
neuromechIntegratedVulnerability( emg, kinematics, force_data, hg = NULL, emg_mapping = NULL, kin_mapping = NULL, weights = c(1, 1, 1), use_proxy = TRUE )neuromechIntegratedVulnerability( emg, kinematics, force_data, hg = NULL, emg_mapping = NULL, kin_mapping = NULL, weights = c(1, 1, 1), use_proxy = TRUE )
emg |
EMG data: SummarizedExperiment, matrix, or vector. |
kinematics |
Kinematic data: matrix (time x segments) or SummarizedExperiment. |
force_data |
Force data: named numeric vector, data.frame, or matrix. |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
kin_mapping |
Optional pre-computed data.frame from |
weights |
Numeric vector of length 3, weights for EMG, kinematics, force components (default: c(1, 1, 1)). |
use_proxy |
Logical, use degree-based proxy for impact deviation (default: TRUE). |
A list with:
Named numeric vector of muscle vulnerability scores
Data.frame sorted by vulnerability (descending)
Pairwise correlations between EMG, kin, force sources
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(4) kin <- matrix(cumsum(rnorm(300)) + sin(seq(0, 6, length.out = 300)), 300, 1) colnames(kin) <- "elbow" set.seed(5) force_data <- matrix(abs(rnorm(300)) + 1, 300, 1); colnames(force_data) <- "grf" result <- neuromechIntegratedVulnerability(emg, kin, force_data) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(4) kin <- matrix(cumsum(rnorm(300)) + sin(seq(0, 6, length.out = 300)), 300, 1) colnames(kin) <- "elbow" set.seed(5) force_data <- matrix(abs(rnorm(300)) + 1, 300, 1); colnames(force_data) <- "grf" result <- neuromechIntegratedVulnerability(emg, kin, force_data) result
Estimates joint torques from EMG activation and anatomical moment arms, computing per-muscle contributions, coactivation indices, and torque balance ratios across joints.
neuromechJointTorque( emg, hg = NULL, emg_mapping = NULL, moment_arm_table = NULL, activation_method = c("rms", "mean_rectified", "peak"), sr = NULL, joints = NULL, n_perm = 999L )neuromechJointTorque( emg, hg = NULL, emg_mapping = NULL, moment_arm_table = NULL, activation_method = c("rms", "mean_rectified", "peak"), sr = NULL, joints = NULL, n_perm = 999L )
emg |
EMG data: SummarizedExperiment, matrix (time x channels), or vector. |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
moment_arm_table |
Optional custom data.frame with columns: muscle_name, joint_name, moment_arm_m, direction. Overrides the built-in lookup. |
activation_method |
Character, method for computing activation: "rms" (default), "mean_rectified", "peak". |
sr |
Optional sampling rate. |
joints |
Optional character vector of joint names to restrict analysis. |
n_perm |
Integer, number of permutations for Mantel test (default: 999). |
An S3 object of class "MSKNeuromechTorque" with:
Data.frame: muscle, joint, activation, moment_arm, direction, torque_contribution
Data.frame: joint, net_torque, agonist_sum, antagonist_sum, coactivation_index, n_muscles
Data.frame: joint, balance_ratio
Mantel test result
Character: "lookup" or "custom"
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- neuromechJointTorque(emg, sr = 1000, n_perm = 19) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- neuromechJointTorque(emg, sr = 1000, n_perm = 19) result
Maps cortical drive (CMC), muscle activation (EMG RMS), and force output to the MSK hypergraph to compute per-muscle drive efficiency and aggregate per community.
neuromechMotorDriveTopography( eeg, emg, force_data, hg = NULL, freq_band = c(15, 35), gamma = 4.3, emg_mapping = NULL, sr = NULL )neuromechMotorDriveTopography( eeg, emg, force_data, hg = NULL, freq_band = c(15, 35), gamma = 4.3, emg_mapping = NULL, sr = NULL )
eeg |
EEG data: SummarizedExperiment, matrix, or vector. |
emg |
EMG data: SummarizedExperiment, matrix, or vector. |
force_data |
Force data: named numeric vector (per muscle), data.frame with columns (muscle, force), or matrix (time x channels for RMS). |
hg |
An MSKHypergraph object (NULL loads default). |
freq_band |
Numeric vector of length 2, CMC frequency band (default: c(15, 35)). |
gamma |
Numeric, resolution parameter for community detection (default: 4.3). |
emg_mapping |
Optional pre-computed data.frame from |
sr |
Optional sampling rate. |
A list with:
Data.frame with muscle, cortical_drive, emg_activation, force, drive_efficiency
Data.frame with community, mean drive/activation/force
Pearson r between cortical drive and force
p-value for the correlation
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") set.seed(5) force_data <- matrix(abs(rnorm(300)) + 1, 300, 1); colnames(force_data) <- "grf" result <- neuromechMotorDriveTopography(eeg, emg, force_data, sr = 1000) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") set.seed(5) force_data <- matrix(abs(rnorm(300)) + 1, 300, 1); colnames(force_data) <- "grf" result <- neuromechMotorDriveTopography(eeg, emg, force_data, sr = 1000) result
Extracts muscle synergies from EMG data via NMF or PCA, maps them to MSK structural communities, and compares synergy-based partitioning with structural communities via z-Rand.
neuromechMuscleSynergy( emg, hg = NULL, emg_mapping = NULL, method = c("nmf", "pca"), n_synergies = 4L, auto_select = FALSE, vaf_threshold = 0.9, max_k = 10L, gamma = 4.3, n_perm = 999L, seed = NULL, sr = NULL )neuromechMuscleSynergy( emg, hg = NULL, emg_mapping = NULL, method = c("nmf", "pca"), n_synergies = 4L, auto_select = FALSE, vaf_threshold = 0.9, max_k = 10L, gamma = 4.3, n_perm = 999L, seed = NULL, sr = NULL )
emg |
EMG data: SummarizedExperiment, matrix (time x channels), or vector. |
hg |
An MSKHypergraph object (NULL loads default). |
emg_mapping |
Optional pre-computed data.frame from |
method |
Character, decomposition method: "nmf" (default) or "pca". |
n_synergies |
Integer, number of synergies to extract (default: 4). |
auto_select |
Logical, automatically select n_synergies via VAF criterion (default: FALSE). |
vaf_threshold |
Numeric, VAF threshold for auto selection (default: 0.90). |
max_k |
Integer, maximum k to try during auto selection (default: 10). |
gamma |
Numeric, resolution parameter for MSK community detection (default: 4.3). |
n_perm |
Integer, permutations for z-Rand (unused, reserved). |
seed |
Optional integer seed for NMF reproducibility. |
sr |
Optional sampling rate. |
An S3 object of class "MSKNeuromechSynergy" with:
Synergy weight matrix (n_muscles x n_synergies)
Activation coefficients (n_synergies x n_timepoints)
Number of synergies extracted
Variance accounted for
VAF curve if auto_select (else NULL)
Decomposition method used
Synergy-to-community enrichment data.frame
Cosine similarity matrix (k x k)
z-Rand comparing synergy vs structural
Reconstruction error
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- neuromechMuscleSynergy(emg, method = "nmf", n_synergies = 2, n_perm = 19, seed = 1, sr = 1000) resultset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") result <- neuromechMuscleSynergy(emg, method = "nmf", n_synergies = 2, n_perm = 19, seed = 1, sr = 1000) result
Orchestrator function that runs all available neuromechanics analyses with shared mapping, skipping analyses when inputs are NULL.
neuromechSummary( eeg = NULL, emg, kinematics = NULL, force_data = NULL, hg = NULL, freq_band = c(15, 35), gamma = 4.3, sr = NULL, n_synergies = 4L, synergy_method = "nmf" )neuromechSummary( eeg = NULL, emg, kinematics = NULL, force_data = NULL, hg = NULL, freq_band = c(15, 35), gamma = 4.3, sr = NULL, n_synergies = 4L, synergy_method = "nmf" )
eeg |
Optional EEG data (NULL to skip CMC and motor drive analyses). |
emg |
EMG data (required). |
kinematics |
Optional kinematic data (NULL to skip EMD and vulnerability). |
force_data |
Optional force data (NULL to skip motor drive and vulnerability). |
hg |
An MSKHypergraph object (NULL loads default). |
freq_band |
Numeric vector of length 2, CMC frequency band (default: c(15, 35)). |
gamma |
Numeric, resolution parameter (default: 4.3). |
sr |
Optional sampling rate. |
n_synergies |
Integer, number of synergies for muscle synergy analysis (default: 4). |
synergy_method |
Character, synergy decomposition method: "nmf" (default) or "pca". |
An S3 object of class "MSKNeuromechSummary" with:
CMC result or NULL
EMD result or NULL
Motor drive result or NULL
Vulnerability result or NULL
Joint torque result or NULL
Muscle synergy result or NULL
Directional coupling result or NULL
Character vector of successfully computed analyses
set.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") set.seed(4) kin <- matrix(cumsum(rnorm(300)) + sin(seq(0, 6, length.out = 300)), 300, 1) colnames(kin) <- "elbow" set.seed(5) force_data <- matrix(abs(rnorm(300)) + 1, 300, 1); colnames(force_data) <- "grf" s <- neuromechSummary(eeg = eeg, emg = emg, kinematics = kin, force_data = force_data, sr = 1000, n_synergies = 2) sset.seed(1) emg <- matrix(abs(rnorm(300 * 4)), 300, 4) colnames(emg) <- c("Biceps Brachii", "Deltoid", "Trapezius", "Triceps Brachii") set.seed(3) eeg <- matrix(rnorm(300 * 3), 300, 3); colnames(eeg) <- c("C3", "Cz", "C4") set.seed(4) kin <- matrix(cumsum(rnorm(300)) + sin(seq(0, 6, length.out = 300)), 300, 1) colnames(kin) <- "elbow" set.seed(5) force_data <- matrix(abs(rnorm(300)) + 1, 300, 1); colnames(force_data) <- "grf" s <- neuromechSummary(eeg = eeg, emg = emg, kinematics = kin, force_data = force_data, sr = 1000, n_synergies = 2) s
Runs network analysis on multiple OpenSim models and performs statistical comparisons of their network properties.
opensimCompareModels(model_paths, names = NULL, gamma = 4.3)opensimCompareModels(model_paths, names = NULL, gamma = 4.3)
model_paths |
Character vector of paths to .osim files. |
names |
Optional character vector of model names. |
gamma |
Resolution parameter for community detection (default: 4.3). |
Requires a native OpenSim installation and an .osim model file,
so the example cannot run offline and is not executed.
An S3 object of class "MSKModelComparison" with:
List of MSKOpenSimAnalysis objects
Data frame of pairwise comparison statistics
## Not run: # Requires a native OpenSim install and an .osim model file. result <- opensimCompareModels(c("modelA.osim", "modelB.osim")) ## End(Not run)## Not run: # Requires a native OpenSim install and an .osim model file. result <- opensimCompareModels(c("modelA.osim", "modelB.osim")) ## End(Not run)
Combines muscle force data from OpenSim simulations with MSK network impact scores to compute force-weighted vulnerability.
opensimForceToImpact( model_path, force_data, hg = NULL, method = c("weighted_simulation", "weighted_scores") )opensimForceToImpact( model_path, force_data, hg = NULL, method = c("weighted_simulation", "weighted_scores") )
model_path |
Character, path to .osim file. |
force_data |
Force data as a data.frame (with muscle name column), named numeric vector, or path to an OpenSim .sto file. |
hg |
An optional pre-built MSKHypergraph (if NULL, built from model). |
method |
Character, weighting method:
|
Requires a native OpenSim installation and an .osim model file,
so the example cannot run offline and is not executed.
A list with:
Force-weighted impact scores
Force values used for weighting
Original (unweighted) impact scores
## Not run: # Requires a native OpenSim install and an .osim model file. result <- opensimForceToImpact("gait2392.osim", force_data = matrix(0, 100, 3)) ## End(Not run)## Not run: # Requires a native OpenSim install and an .osim model file. result <- opensimForceToImpact("gait2392.osim", force_data = matrix(0, 100, 3)) ## End(Not run)
Builds a hypergraph from an OpenSim model and runs complete network analysis including metrics, community detection, and optionally impact scoring.
opensimNetworkAnalysis( model = NULL, model_path = NULL, gamma = 4.3, run_simulation = TRUE )opensimNetworkAnalysis( model = NULL, model_path = NULL, gamma = 4.3, run_simulation = TRUE )
model |
An optional PhysioOpenSimModel object. |
model_path |
Character, path to .osim file. |
gamma |
Resolution parameter for community detection (default: 4.3). |
run_simulation |
Logical, whether to run full impact simulation (default: TRUE; set to FALSE for faster analysis). |
Requires a native OpenSim installation and an .osim model file,
so the example cannot run offline and is not executed.
An S3 object of class "MSKOpenSimAnalysis" with:
The MSKHypergraph object
Network metrics from mskNetworkMetrics
Community detection results
Impact scores (if run_simulation = TRUE)
## Not run: # Requires a native OpenSim install and an .osim model file. result <- opensimNetworkAnalysis(model_path = "gait2392.osim") ## End(Not run)## Not run: # Requires a native OpenSim install and an .osim model file. result <- opensimNetworkAnalysis(model_path = "gait2392.osim") ## End(Not run)
Parses an OpenSim .osim file to extract muscle-bone attachment points and constructs a subject-specific MSKHypergraph.
opensimToMSKHypergraph(model = NULL, model_path = NULL)opensimToMSKHypergraph(model = NULL, model_path = NULL)
model |
An optional PhysioOpenSimModel object. If provided, the model's path is used and results are cross-validated. |
model_path |
Character, path to a .osim file. Required if |
Requires a native OpenSim installation and an .osim model file,
so the example cannot run offline and is not executed.
An MSKHypergraph object with subject-specific anatomy.
Requires the xml2 package.
## Not run: hg <- opensimToMSKHypergraph(model_path = "gait2392.osim") print(hg) ## End(Not run)## Not run: hg <- opensimToMSKHypergraph(model_path = "gait2392.osim") print(hg) ## End(Not run)
Visualizes detected communities with size information.
plotCommunityStructure(community, ...)plotCommunityStructure(community, ...)
community |
Result from mskCommunityDetect or mskConsensusPartition. |
... |
Additional arguments. |
Invisible NULL.
cm <- mskCommunityDetect(MSKHypergraph(), gamma = 4.3) plotCommunityStructure(cm)cm <- mskCommunityDetect(MSKHypergraph(), gamma = 4.3) plotCommunityStructure(cm)
Plots the degree distribution of the hypergraph, comparing real and null model distributions. Reproduces Fig 2e.
plotDegreeDistribution( hg = NULL, type = c("muscle", "bone"), null_dist = NULL, log_scale = TRUE, ... )plotDegreeDistribution( hg = NULL, type = c("muscle", "bone"), null_dist = NULL, log_scale = TRUE, ... )
hg |
An MSKHypergraph object. |
type |
"muscle" or "bone" distribution. |
null_dist |
Optional numeric vector of null model degree counts. |
log_scale |
Logical, use log scale (default: TRUE for heavy-tail). |
... |
Additional arguments passed to barplot. |
Invisible NULL.
hg <- MSKHypergraph() plotDegreeDistribution(hg, type = "muscle")hg <- MSKHypergraph() plotDegreeDistribution(hg, type = "muscle")
Recreates Fig 4b: deviation ratio vs homunculus area. Target: F(1,19) = 21.3, R^2 = 0.52
plotHomunculus(homunculus_data = NULL, ...)plotHomunculus(homunculus_data = NULL, ...)
homunculus_data |
Optional data.frame. If NULL, uses built-in data. |
... |
Additional arguments. |
Invisible the regression result.
plotHomunculus()plotHomunculus()
Recreates Fig 3a: scatter plot of impact scores by hyperedge degree. Target: R^2 = 0.45 for the relationship.
plotImpactVsDegree(impact_scores, hg = NULL, show_regression = TRUE, ...)plotImpactVsDegree(impact_scores, hg = NULL, show_regression = TRUE, ...)
impact_scores |
Named numeric vector of impact scores. |
hg |
An MSKHypergraph object. |
show_regression |
Logical, overlay regression line. |
... |
Additional arguments. |
Invisible the regression result.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) scores <- mskImpactScoreAll(mskSimulate(hg, n_steps = 50), verbose = FALSE) plotImpactVsDegree(scores, hg, show_regression = FALSE)C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) scores <- mskImpactScoreAll(mskSimulate(hg, n_steps = 50), verbose = FALSE) plotImpactVsDegree(scores, hg, show_regression = FALSE)
Recreates Fig 3b: weighted regression of impact deviation vs clinical recovery time. Target: R^2 = 0.757.
plotImpactVsRecovery(recovery_data = NULL, ...)plotImpactVsRecovery(recovery_data = NULL, ...)
recovery_data |
Optional data.frame. If NULL, uses built-in data. |
... |
Additional arguments. |
Invisible the regression result.
plotImpactVsRecovery()plotImpactVsRecovery()
Visualizes the musculoskeletal network with optional community coloring.
plotMSKNetwork( hg = NULL, type = c("muscle", "bone"), membership = NULL, layout = "fr", vertex_size = 3, ... )plotMSKNetwork( hg = NULL, type = c("muscle", "bone"), membership = NULL, layout = "fr", vertex_size = 3, ... )
hg |
An MSKHypergraph object. |
type |
"bone" or "muscle" projection to plot. |
membership |
Optional community membership vector for coloring. |
layout |
Character, layout algorithm: "fr" (Fruchterman-Reingold), "circle", or "auto". |
vertex_size |
Numeric, base vertex size. |
... |
Additional arguments passed to plot. |
Invisible NULL. Produces a plot.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) plotMSKNetwork(MSKHypergraph(C), type = "bone")C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) plotMSKNetwork(MSKHypergraph(C), type = "bone")
Creates the bone-bone weighted adjacency matrix A = t(C) %*% C where A[i,j] counts the number of muscles shared between bones i and j. This is the one-mode projection onto the bone (vertex) space.
projectBoneGraph(hg)projectBoneGraph(hg)
hg |
An MSKHypergraph object, or a sparse incidence matrix. |
A symmetric sparse Matrix (n_bones x n_bones)
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) as.matrix(projectBoneGraph(hg))C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) as.matrix(projectBoneGraph(hg))
Creates the muscle-muscle weighted adjacency matrix B = C %*% t(C) where B[i,j] counts the number of bones shared between muscles i and j. This is the one-mode projection onto the muscle (hyperedge) space.
projectMuscleGraph(hg)projectMuscleGraph(hg)
hg |
An MSKHypergraph object, or a sparse incidence matrix. |
A symmetric sparse Matrix (n_muscles x n_muscles)
Murphy AC et al. (2018) PLOS Biology.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) as.matrix(projectMuscleGraph(hg))C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) hg <- MSKHypergraph(C) as.matrix(projectMuscleGraph(hg))
Normalised tendon force as a function of normalised tendon length
(tendon length / tendon slack length): zero below slack, rising steeply to
one maximum-isometric-force at the reference strain.
tendonForceLength(norm_tendon_length, e0t = 0.04, kt = 3)tendonForceLength(norm_tendon_length, e0t = 0.04, kt = 3)
norm_tendon_length |
Tendon length divided by tendon slack length. |
e0t |
Tendon strain at maximum isometric force (default 0.04). |
kt |
Exponential stiffness shape (default 3). |
Normalised tendon force (>= 0).
tendonForceLength(c(1.0, 1.02, 1.04))tendonForceLength(c(1.0, 1.02, 1.04))
Returns the degree of each vertex (bone), i.e., the number of muscles attached to each bone.
vertexDegree(hg)vertexDegree(hg)
hg |
An MSKHypergraph object or incidence matrix. |
Named numeric vector of bone degrees.
C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) vertexDegree(MSKHypergraph(C))C <- matrix(c(1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 4, dimnames = list(paste0("b", 1:4), paste0("m", 1:4))) vertexDegree(MSKHypergraph(C))