| Title: | Digital-Twin Simulation and Data Assimilation for Physiological and Movement Systems |
|---|---|
| Description: | A digital-twin engine that couples mechanistic simulation of movement and physiology with realistic sensor emulation and data assimilation, so a subject-specific model can be personalised from real measurements and then used for in-silico experiments. It ships a parameterised limb model (forward dynamics) and mechanistic physiological and network generators -- a Jansen-Rit neural-mass EEG column with an optional 1/f aperiodic background, a McSharry limit-cycle electrocardiogram, a Windkessel arterial blood-pressure model with a closed-loop baroreflex, a respiration waveform and respiratory-sinus-arrhythmia heart-rate-variability generator, and Kuramoto coupled-oscillator and Matsuoka central-pattern-generator rhythms. A forward sensor-emulation layer turns clean signals into realistic recordings (inertial-measurement-unit, optical-marker occlusion, surface-electromyography electrode and electrocardiogram lead), with richer chains for multi-channel electromyography crosstalk and amplifier characteristics, multi-lead electrocardiogram and analogue-to-digital conversion. A data-assimilation layer personalises the twin to recorded data with a family of state-space estimators (unscented, extended, particle and ensemble Kalman filters), likelihood-free Bayesian calibration (rejection and sequential Monte-Carlo approximate Bayesian computation), Markov-chain Monte- Carlo, a Gaussian-process surrogate and GP-accelerated (Bayesian-optimisation) calibration -- including personalising the closed-loop cardiovascular- respiratory twin to measured heart-rate and blood-pressure variability -- as well as particle Markov-chain Monte-Carlo and particle Gibbs (with ancestor sampling) for whole-waveform Bayesian fitting -- including a whole-record fit of the closed-loop cardiovascular model that identifies the baroreflex and respiratory gains a summary fit cannot, with either a resonant Mayer-oscillator or a broadband first-order-autoregressive low-frequency component for records without a distinct Mayer peak -- and Fisher-information optimal experimental design for choosing the most informative experimental condition. Applications include single- and population-level in-silico intervention (perturb parameters and predict the change, with effect size, responder rate and statistical power), multi-modal coupling, and labelled synthetic training-data generation with sensor domain randomisation, all backed by a simulation validity-and-verification harness (identifiability, global sensitivity and validation against real recordings). A clinical decision loop ties these together as research tooling: a personalised twin is certified by a validation gate (goodness of fit and parameter identifiability) before it may advise, candidate interventions are scored with their predictive uncertainty and ranked on multiple criteria, mechanistic predictions are reconciled with external evidence rather than overriding it, and a follow-up measurement is checked against the predicted interval. A closed-loop cardiovascular-respiratory model couples a delayed baroreflex with respiration to produce emergent Mayer waves and respiratory sinus arrhythmia, and an optional bridge holds the twin's predictions to account against a real single-case (SCED/MCID/ICF) clinical evaluation. |
| Authors: | Yusuke Matsui [aut, cre] |
| Maintainer: | Yusuke Matsui <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.2 |
| Built: | 2026-10-03 13:43:09 UTC |
| Source: | https://github.com/x-biosignal/PhysioTwin |
Calibrates a simulator's parameters to observed data without a likelihood:
parameters are drawn from prior, passed through simulate to produce summary
statistics, and the draws whose summaries fall closest to observed (the
smallest-distance accept fraction) are retained as posterior samples.
Distances are Euclidean on summaries scaled by their simulated SD, so
statistics on different scales contribute comparably. With regression = TRUE
the accepted draws are corrected by a distance-weighted local-linear regression
onto the observed summaries.
abcCalibration( simulate, observed, prior, n = 5000, accept = 0.02, regression = TRUE, seed = NULL )abcCalibration( simulate, observed, prior, n = 5000, accept = 0.02, regression = TRUE, seed = NULL )
simulate |
|
observed |
Observed summary statistics (same length as |
prior |
|
n |
Number of prior simulations. |
accept |
Fraction of draws to retain (the ABC tolerance as a quantile). |
regression |
If |
seed |
Optional RNG seed. |
an abc_posterior: samples (accepted, adjusted parameter draws),
raw (accepted draws before adjustment), mean, sd, quantiles (2.5/50/
97.5%), threshold (accepted distance) and accept_rate.
Beaumont MA, Zhang W, Balding DJ (2002) Genetics 162:2025-2035; Pritchard JK et al. (1999).
personalizeTwin(), profileLikelihood(), gpEmulator()
# recover the mean/SD of a Gaussian from its sample mean and SD set.seed(1); obs <- c(mean = 2, sd = 1.5) sim <- function(th) { x <- rnorm(200, th[1], th[2]); c(mean(x), sd(x)) } pri <- function() c(m = runif(1, -3, 7), s = runif(1, 0.2, 4)) post <- abcCalibration(sim, obs, pri, n = 2000, accept = 0.02, seed = 1) post$mean# recover the mean/SD of a Gaussian from its sample mean and SD set.seed(1); obs <- c(mean = 2, sd = 1.5) sim <- function(th) { x <- rnorm(200, th[1], th[2]); c(mean(x), sd(x)) } pri <- function() c(m = runif(1, -3, 7), s = runif(1, 0.2, 4)) post <- abcCalibration(sim, obs, pri, n = 2000, accept = 0.02, seed = 1) post$mean
A sequential Monte-Carlo approximate Bayesian computation sampler. A population of parameter particles is carried through a schedule of shrinking tolerances: at each level a particle is drawn from the previous population, perturbed by a Gaussian kernel (covariance twice the population's), simulated, and accepted if its summary distance is below the current tolerance; importance weights correct for the kernel. The tolerance for each level is an adaptive quantile of the previous population's distances. This concentrates on the posterior with far fewer simulations than one-shot rejection ABC.
abcSMC( simulate, observed, prior, prior_density = NULL, n_particles = 200, n_populations = 5, quantile = 0.5, max_tries = 2000L, seed = NULL )abcSMC( simulate, observed, prior, prior_density = NULL, n_particles = 200, n_populations = 5, quantile = 0.5, max_tries = 2000L, seed = NULL )
simulate |
|
observed |
Observed summary statistics. |
prior |
|
prior_density |
Optional |
n_particles |
Particles per population. |
n_populations |
Number of tolerance levels. |
quantile |
Tolerance schedule: each level's tolerance is this quantile of the previous population's distances. |
max_tries |
Per-particle acceptance attempts before giving up (guards a too-tight tolerance). |
seed |
Optional RNG seed. |
an abc_smc: weighted posterior samples, weights, the tolerance
schedule, and the weighted mean/sd.
Toni T et al. (2009) J R Soc Interface 6:187-202; Beaumont MA et al. (2009) Biometrika 96:983-990.
abcCalibration(), metropolis()
set.seed(1); obs <- c(mean = 2, sd = 1.5) sim <- function(th) { x <- rnorm(200, th[1], th[2]); c(mean(x), sd(x)) } pri <- function() c(m = runif(1, -3, 7), s = runif(1, 0.2, 4)) post <- abcSMC(sim, obs, pri, n_particles = 100, n_populations = 4, seed = 1) post$meanset.seed(1); obs <- c(mean = 2, sd = 1.5) sim <- function(th) { x <- rnorm(200, th[1], th[2]); c(mean(x), sd(x)) } pri <- function() c(m = runif(1, -3, 7), s = runif(1, 0.2, 4)) post <- abcSMC(sim, obs, pri, n_particles = 100, n_populations = 4, seed = 1) post$mean
Rounds a signal to bits resolution over range, as an analogue-to-digital
converter does. Optional (subtractive) dither randomises the quantisation error,
decorrelating it from the signal.
adcQuantize(signal, bits = 12, range = NULL, dither = FALSE)adcQuantize(signal, bits = 12, range = NULL, dither = FALSE)
signal |
Numeric signal. |
bits |
Bit depth (e.g. 8, 12, 16). |
range |
|
dither |
If |
the quantised signal (clipped to range).
adcSample(), quantizationNoise()
q <- adcQuantize(sin(seq(0, 6, length.out = 500)), bits = 8) length(unique(round(q, 6))) <= 256q <- adcQuantize(sin(seq(0, 6, length.out = 500)), bits = 8) length(unique(round(q, 6))) <= 256
Resamples signal from fs_in to fs_out by linear interpolation, optionally
perturbing each output sample time by Gaussian timing jitter – the aperture
jitter of a real converter's clock.
adcSample(signal, fs_in, fs_out, jitter = 0)adcSample(signal, fs_in, fs_out, jitter = 0)
signal |
Numeric signal sampled at |
fs_in, fs_out
|
Input and output sampling rates (Hz). |
jitter |
Sampling-time jitter SD, in input samples. |
a list: signal (resampled), fs (= fs_out), time.
down <- adcSample(sin(seq(0, 12, length.out = 2000)), fs_in = 1000, fs_out = 250) down$fsdown <- adcSample(sin(seq(0, 12, length.out = 2000)), fs_in = 1000, fs_out = 250) down$fs
Generates zero-mean noise whose power spectrum falls as 1 / f^exponent – the
aperiodic ("1/f") background that underlies real electrophysiological recordings
and on which oscillatory rhythms sit. It is made by shaping white noise in the
frequency domain. Added to a mechanistic oscillator (e.g. jansenRit()), it
turns an over-concentrated single-rhythm spectrum into a realistic broadband one.
aperiodicNoise(n, exponent = 1.5, sd = 1, seed = NULL)aperiodicNoise(n, exponent = 1.5, sd = 1, seed = NULL)
n |
Number of samples. |
exponent |
Spectral slope: power scales as |
sd |
Target standard deviation of the output. |
seed |
Optional RNG seed. |
a length-n numeric vector.
He BJ (2014) Scale-free brain activity. Trends Cogn Sci 18:480-487; Donoghue T et al. (2020) Nat Neurosci 23:1655-1665 (aperiodic/periodic).
x <- aperiodicNoise(1024, exponent = 1.5, seed = 1) round(sd(x), 3)x <- aperiodicNoise(1024, exponent = 1.5, seed = 1) round(sd(x), 3)
Convenience helper: the fraction of spectral power of a jansenRit EEG (after
dropping the transient) that falls in each conventional band.
bandPower(x)bandPower(x)
x |
a |
a named numeric vector of relative power in delta, theta, alpha,
beta and gamma.
bandPower(jansenRit(duration = 6, seed = 1))["alpha"] # alpha-dominantbandPower(jansenRit(duration = 6, seed = 1))["alpha"] # alpha-dominant
Closes the short-term baroreflex around the windkessel() arterial model:
each beat the mean arterial pressure is compared to a set-point and the next
heart period (RR interval) is set to RR0 + gain * (MAP - setpoint) – a
higher pressure lengthens RR (slows the heart), lowering cardiac output and
buffering the pressure (negative feedback; gain is the baroreflex
sensitivity, s/mmHg). A step change in peripheral resistance at perturb_time
(e.g. vasodilation) challenges the loop; with the reflex engaged the mean
pressure is actively regulated back toward the set-point.
baroreflex( duration = 60, setpoint = 92, gain = 0.02, hr0 = 70, stroke_volume = 80, R = 1, C = 1.9, Zc = 0.05, systolic_frac = 0.3, perturb_time = 30, perturb_R = 0.7, reflex = TRUE, sfp = 250 )baroreflex( duration = 60, setpoint = 92, gain = 0.02, hr0 = 70, stroke_volume = 80, R = 1, C = 1.9, Zc = 0.05, systolic_frac = 0.3, perturb_time = 30, perturb_R = 0.7, reflex = TRUE, sfp = 250 )
duration |
Recording length (s). |
setpoint |
Regulated mean-arterial-pressure set-point (mmHg). |
gain |
Baroreflex sensitivity, RR change per unit pressure error (s/mmHg). |
hr0 |
Baseline heart rate (bpm). |
stroke_volume, R, C, Zc, systolic_frac
|
Windkessel parameters (see |
perturb_time |
Time of the resistance step (s); |
perturb_R |
Multiplier applied to |
reflex |
If |
sfp |
Output sampling rate (Hz). |
a baroreflex object: continuous time/pressure; per-beat
beat_time, rr, hr, beat_map; the setpoint; and regulated (the
post-perturbation steady-state absolute MAP error).
Ursino M (1998) Am J Physiol 275:H1733-H1747. deBoer RW et al. (1987) Am J Physiol 253:H680-H689.
br <- baroreflex(duration = 60, perturb_time = 30, perturb_R = 0.7) br$regulated # small residual MAP error -- pressure is bufferedbr <- baroreflex(duration = 60, perturb_time = 30, perturb_R = 0.7) br$regulated # small residual MAP error -- pressure is buffered
A beat-to-beat model in which arterial pressure, heart period and peripheral resistance evolve under a delayed, two-arm baroreflex, coupled to respiration. A Windkessel map sets each beat's systolic and diastolic pressure; the baroreflex adjusts the next heart period (fast vagal + delayed sympathetic) and the peripheral resistance (delayed sympathetic) from the sensed systolic pressure; respiration modulates stroke volume (a mechanical blood-pressure effect) and vagal tone (respiratory sinus arrhythmia). The delayed resistance/heart-period feedback, excited by small noise, produces an EMERGENT ~0.1 Hz Mayer wave in pressure and heart rate; respiration adds the respiratory-frequency oscillations. None of these rhythms is imposed.
cardioRespiratory( duration = 300, hr0 = 70, resp_rate = 12, S0 = 120, D0 = 80, R0 = 1, g_vagal = 0.01, g_symp = 0.007, symp_delay = 2L, g_resist = 0.01, resist_delay = 4L, rsa_amp = 0.035, resp_bp_amp = 2.5, starling = 0.1, noise = 0.02, respiration = TRUE, seed = NULL )cardioRespiratory( duration = 300, hr0 = 70, resp_rate = 12, S0 = 120, D0 = 80, R0 = 1, g_vagal = 0.01, g_symp = 0.007, symp_delay = 2L, g_resist = 0.01, resist_delay = 4L, rsa_amp = 0.035, resp_bp_amp = 2.5, starling = 0.1, noise = 0.02, respiration = TRUE, seed = NULL )
duration |
Recording length (s). |
hr0 |
Baseline heart rate (bpm). |
resp_rate |
Breathing rate (breaths/min). |
S0, D0
|
Target systolic/diastolic pressure (mmHg). |
R0 |
Baseline peripheral resistance (relative). |
g_vagal, g_symp
|
Baroreflex heart-period gains (s/mmHg): fast vagal (current beat) and delayed sympathetic. |
symp_delay, resist_delay
|
Sympathetic heart-period and resistance baroreflex delays (beats). |
g_resist |
Baroreflex resistance gain (per mmHg). |
rsa_amp |
Respiratory vagal modulation of heart period (s). |
resp_bp_amp |
Respiratory mechanical modulation of stroke volume (mmHg). |
starling |
Frank-Starling sensitivity of stroke volume to the preceding heart period (fractional). |
noise |
SD of the beat-to-beat neural noise (fraction) that excites the loop. |
respiration |
If |
seed |
Optional RNG seed. |
a cardio_respiratory object: per-beat time, rr, hr, systolic,
diastolic, resistance; the HRV and systolic-pressure LF/HF band summaries
hrv and bpv (each with lf_power/hf_power/lf_peak/hf_peak); and the
baroreflex sensitivity brs (ms/mmHg) recovered from the beat series.
de Boer RW, Karemaker JM, Strackee J (1987) Am J Physiol 253:H680-H689; Julien C (2006) Cardiovasc Res 70:12-21 (Mayer waves).
windkessel(), baroreflex(), rsaTachogram()
cr <- cardioRespiratory(duration = 240, resp_rate = 12, seed = 1) cr$bpv$lf_peak # ~0.1 Hz Mayer wave in blood pressure cr$hrv$hf_peak # respiratory-frequency RSA in heart ratecr <- cardioRespiratory(duration = 240, resp_rate = 12, seed = 1) cr$bpv$lf_peak # ~0.1 Hz Mayer wave in blood pressure cr$hrv$hf_peak # respiratory-frequency RSA in heart rate
A linear-Gaussian state-space representation of the RR-interval series: a
damped second-order Mayer oscillator (the delayed-baroreflex resonance) drives
the heart period through a gain g, and respiration adds a high-frequency
(respiratory sinus arrhythmia) term of amplitude A_rsa. Suitable as the
build_model argument of particleMCMC() / particleGibbs().
cardioStateSpace( g = 8, rho = 0.9, omega_m = 0.53, A_rsa = 20, resp, rr0 = 850, q_mayer = 1, r_obs = 9, lf_model = c("oscillator", "ar1"), phi = 0.9 )cardioStateSpace( g = 8, rho = 0.9, omega_m = 0.53, A_rsa = 20, resp, rr0 = 850, q_mayer = 1, r_obs = 9, lf_model = c("oscillator", "ar1"), phi = 0.9 )
g |
Low-frequency gain onto the RR interval (ms per unit state). |
rho |
Mayer-oscillator damping (0..1; nearer 1 = sharper low-frequency peak;
used when |
omega_m |
Mayer angular frequency (radians per beat; |
A_rsa |
Respiratory (RSA) amplitude on the RR interval (ms). |
resp |
Respiratory regressor sampled at the beats (length = number of beats). |
rr0 |
Mean RR interval (ms). |
q_mayer, r_obs
|
Low-frequency process-noise and RR observation-noise variances. |
lf_model |
Low-frequency model: |
phi |
First-order autoregressive coefficient (0..1; used when
|
a state-space model list f, h (function(x, k)), Q, R, x0, P0.
personalizeCardioWaveform(), particleMCMC(), cardioRespiratory()
resp <- sin(2 * pi * 0.25 * cumsum(rep(0.85, 200))) m <- cardioStateSpace(g = 8, rho = 0.9, omega_m = 0.53, A_rsa = 20, resp = resp) m$h(c(1, 0), 1)resp <- sin(2 * pi * 0.25 * cumsum(rep(0.85, 200))) m <- cardioStateSpace(g = 8, rho = 0.9, omega_m = 0.53, A_rsa = 20, resp = resp) m$h(c(1, 0), 1)
The heart-rate/blood-pressure variability summaries used to personalise the
closed-loop twin: baroreflex sensitivity (brs, ms/mmHg), the Mayer-wave
low-frequency power of blood pressure (bp_lf) and the respiratory
high-frequency power of heart rate (hr_hf).
cardioSummary(cr)cardioSummary(cr)
cr |
a |
a named numeric vector c(brs, bp_lf, hr_hf).
cardioRespiratory(), personalizeCardio()
cardioSummary(cardioRespiratory(duration = 120, seed = 1))cardioSummary(cardioRespiratory(duration = 120, seed = 1))
Personalises a twin to a patient's IMU data, gates it with twinReadiness(),
and – only reporting through the gate's confidence – scores and ranks the
candidate interventions. A one-call demonstration of the loop; the returned
object holds the fit, the readiness verdict and the ranked recommendations.
clinicalDecision( imu, prior, estimate = "strength", candidates = defaultInterventions(), outcome = "rom", method = c("optim", "ukf"), dt = 0.005, n_mc = 200, weights = c(efficacy = 1, certainty = 1, risk = 1, burden = 0.5), seed = NULL )clinicalDecision( imu, prior, estimate = "strength", candidates = defaultInterventions(), outcome = "rom", method = c("optim", "ukf"), dt = 0.005, n_mc = 200, weights = c(efficacy = 1, certainty = 1, risk = 1, burden = 0.5), seed = NULL )
imu |
Patient IMU data (as from |
prior |
A |
estimate |
Parameters to personalise. |
candidates |
Candidate interventions (default |
outcome |
Outcome to optimise. |
method |
Personalisation method ( |
dt, n_mc, weights, seed
|
As in |
a clinical_decision: fit, readiness, and scores.
personalizeTwin(), twinReadiness(), scoreInterventions()
truth <- limbTwin(strength = 6, damping = 0.6, stiffness = 3) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) dec <- clinicalDecision(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01, n_mc = 40, seed = 1) dec$scores$ranking$intervention[1]truth <- limbTwin(strength = 6, damping = 0.6, stiffness = 3) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) dec <- clinicalDecision(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01, n_mc = 40, seed = 1) dec$scores$ranking$intervention[1]
Simulates a Matsuoka two-neuron half-centre oscillator: two mutually-inhibiting neurons with fatigue/adaptation produce a self-sustained anti-phase rhythm – the mechanistic building block of locomotor central pattern generators. The output is the difference of the two rectified activities.
cpgMatsuoka( tau = 0.25, tau2 = 0.5, beta = 2.5, gamma = 2.5, u = 1, duration = 20, dt = 0.005 )cpgMatsuoka( tau = 0.25, tau2 = 0.5, beta = 2.5, gamma = 2.5, u = 1, duration = 20, dt = 0.005 )
tau, tau2
|
Rise and adaptation time constants (s). |
beta |
Adaptation (self-inhibition) strength. |
gamma |
Mutual-inhibition strength. |
u |
Tonic drive. |
duration, dt
|
Simulation length and step (s). |
a cpg object: time, output (rhythmic drive y1 - y2), y1, y2
(the two neuron outputs), and frequency (Hz, estimated from the output).
Matsuoka K (1985) Biol Cybern 52:367-376.
cpg <- cpgMatsuoka(duration = 20) cpg$frequency # a sustained locomotor-like rhythmcpg <- cpgMatsuoka(duration = 20) cpg$frequency # a sustained locomotor-like rhythm
A symmetric, energy-normalised mixing matrix: each channel keeps most of its own
signal and picks up crosstalk of an immediate neighbour, crosstalk^2 of the
next, and so on (rows sum to one). crosstalk = 0 gives the identity.
crosstalkMatrix(channels, crosstalk = 0.15)crosstalkMatrix(channels, crosstalk = 0.15)
channels |
Number of channels. |
crosstalk |
Per-neighbour crosstalk fraction (0..1). |
a channels x channels mixing matrix.
round(crosstalkMatrix(3, 0.2), 2)round(crosstalkMatrix(3, 0.2), 2)
A small library of candidate interventions for the single-joint movement twin: strengthening, tone reduction, stretching, an orthosis, and a no-intervention control.
defaultInterventions()defaultInterventions()
a list of twin_intervention objects.
intervention(), scoreInterventions()
names(defaultInterventions())names(defaultInterventions())
High-passes the signal (subtracting a ~1 s moving-average trend) and returns the indices of R waves – local maxima above a fraction of the signal range, separated by a physiological refractory gap.
detectRpeaks(x, fs, refractory = 0.3, thresh = 0.5)detectRpeaks(x, fs, refractory = 0.3, thresh = 0.5)
x |
An ECG signal (numeric). |
fs |
Sampling rate (Hz). |
refractory |
Minimum spacing between R-peaks (s). |
thresh |
Peak threshold as a fraction of the (de-trended) signal range. |
the integer indices of detected R-peaks.
e <- ecgsyn(duration = 10, seed = 1) length(detectRpeaks(e$ecg, e$sfecg))e <- ecgsyn(duration = 10, seed = 1) length(detectRpeaks(e$ecg, e$sfecg))
Bundles the ECG recording chain's artefacts: baseline wander, powerline interference, muscle noise, motion artefact, gain and offset.
ecgLead( gain = 1, noise = 0.01, powerline_hz = 50, powerline_amp = 0.05, baseline_amp = 0.15, baseline_hz = 0.25, motion_amp = 0, offset = 0, fs = 250 )ecgLead( gain = 1, noise = 0.01, powerline_hz = 50, powerline_amp = 0.05, baseline_amp = 0.15, baseline_hz = 0.25, motion_amp = 0, offset = 0, fs = 250 )
gain |
Electrode/amplifier gain. |
noise |
High-frequency muscle-noise SD (mV). |
powerline_hz |
Mains frequency (Hz). |
powerline_amp |
Powerline interference amplitude (mV). |
baseline_amp, baseline_hz
|
Baseline-wander amplitude (mV) and frequency (Hz, typically respiratory ~0.25 Hz). |
motion_amp |
Motion-artefact amplitude (mV; a slow random walk). |
offset |
Electrode DC offset (mV). |
fs |
Sampling rate (Hz). |
an ecg_lead object.
ecgMeasure(), detectRpeaks(), ecgsyn()
ecgLead(baseline_amp = 0.2, powerline_hz = 60)ecgLead(baseline_amp = 0.2, powerline_hz = 60)
Bundles a set of named ECG leads derived from one cardiac source. Each lead has a projection gain and a morphology (derivative-mixing) weight – so leads differ in shape, not only amplitude – plus baseline wander, shared mains interference and muscle noise.
ecgLeadSet( leads = c("II", "V2", "V5"), gains = NULL, baseline_amp = 0.15, baseline_hz = 0.25, powerline_hz = 50, powerline_amp = 0.05, noise = 0.01, fs = 250 )ecgLeadSet( leads = c("II", "V2", "V5"), gains = NULL, baseline_amp = 0.15, baseline_hz = 0.25, powerline_hz = 50, powerline_amp = 0.05, noise = 0.01, fs = 250 )
leads |
Character vector of lead names (e.g. |
gains |
Optional per-lead projection gains; defaults to standard-ish values. |
baseline_amp, baseline_hz
|
Baseline-wander amplitude (mV) and frequency (Hz). |
powerline_hz, powerline_amp
|
Mains frequency (Hz) and amplitude (mV). |
noise |
Muscle-noise SD (mV). |
fs |
Sampling rate (Hz). |
an ecg_lead_set object.
ecgMeasureLeads(), ecgLead(), ecgsyn()
ecgLeadSet(leads = c("II", "V2", "V5"))ecgLeadSet(leads = c("II", "V2", "V5"))
Adds the electrode chain's artefacts – baseline wander, powerline interference, muscle noise, optional motion artefact, gain and offset – to a clean ECG.
ecgMeasure(ecg, sensor = ecgLead(), seed = NULL)ecgMeasure(ecg, sensor = ecgLead(), seed = NULL)
ecg |
Clean ECG signal (mV), e.g. from |
sensor |
An |
seed |
Optional RNG seed. |
an ecg_measurement: signal (recorded lead, mV), clean, time,
fs, sensor.
e <- ecgsyn(duration = 10, seed = 1) rec <- ecgMeasure(e$ecg, ecgLead(baseline_amp = 0.2, fs = e$sfecg), seed = 1) sd(rec$signal)e <- ecgsyn(duration = 10, seed = 1) rec <- ecgMeasure(e$ecg, ecgLead(baseline_amp = 0.2, fs = e$sfecg), seed = 1) sd(rec$signal)
Projects one clean ECG onto a set of leads: each lead is a gain- and morphology-weighted mix of the source and its derivative, plus independent baseline wander, shared mains interference and independent muscle noise.
ecgMeasureLeads(ecg, sensor = ecgLeadSet(), seed = NULL)ecgMeasureLeads(ecg, sensor = ecgLeadSet(), seed = NULL)
ecg |
Clean ECG signal (mV), e.g. from |
sensor |
An |
seed |
Optional RNG seed. |
an ecg_multilead: signals (n_time x n_leads, columns named by lead),
leads, time, fs, sensor.
e <- ecgsyn(duration = 10, seed = 1) ml <- ecgMeasureLeads(e$ecg, ecgLeadSet(fs = e$sfecg), seed = 1) dim(ml$signals)e <- ecgsyn(duration = 10, seed = 1) ml <- ecgMeasureLeads(e$ecg, ecgLeadSet(fs = e$sfecg), seed = 1) dim(ml$signals)
Generates a realistic ECG from the ECGSYN 3-state limit-cycle model: P/Q/R/S/T waves arise as the trajectory passes Gaussian events, and heart-rate variability is injected via a bimodal RR tachogram. The output has physiological morphology AND HRV.
ecgsyn( duration = 10, hr = 60, sfecg = 256, hrv_sd = 0.04, lfhfratio = 0.5, anoise = 0, scale_mv = c(-0.4, 1.2), theta = c(-60, -15, 0, 15, 90) * pi/180, a = c(1.2, -5, 30, -7.5, 0.75), b = c(0.25, 0.1, 0.1, 0.1, 0.4), seed = NULL )ecgsyn( duration = 10, hr = 60, sfecg = 256, hrv_sd = 0.04, lfhfratio = 0.5, anoise = 0, scale_mv = c(-0.4, 1.2), theta = c(-60, -15, 0, 15, 90) * pi/180, a = c(1.2, -5, 30, -7.5, 0.75), b = c(0.25, 0.1, 0.1, 0.1, 0.4), seed = NULL )
duration |
Length (s). |
hr |
Mean heart rate (bpm). |
sfecg |
Sampling frequency (Hz). |
hrv_sd |
RR-interval standard deviation (s) – the amount of HRV (0 = fixed rate). |
lfhfratio |
LF/HF power ratio of the RR spectrum (autonomic balance). |
anoise |
Additive measurement-noise SD (mV). |
scale_mv |
Length-2 physiological output range |
theta, a, b
|
P,Q,R,S,T event angles (rad), amplitudes and widths (McSharry defaults). |
seed |
Optional RNG seed. |
an ecgsyn object: time, ecg, the rr intervals, beats (R-peak
times), sfecg.
McSharry PE, Clifford GD, Tarassenko L, Smith LA (2003) IEEE Trans Biomed Eng 50:289-294.
jansenRit(), windkessel(), rsaTachogram(), kuramoto(), validateTwin()
e <- ecgsyn(duration = 6, hr = 60, sfecg = 256, seed = 1) range(e$ecg) # P/Q/R/S/T morphology length(e$beats) # ~ 6 beatse <- ecgsyn(duration = 6, hr = 60, sfecg = 256, seed = 1) range(e$ecg) # P/Q/R/S/T morphology length(e$beats) # ~ 6 beats
Bundles a multi-channel surface-EMG measurement chain: inter-electrode crosstalk, amplifier gain, a common-mode-rejected mains interference, a DC-blocking high-pass, a passband and thermal noise.
emgArray( channels = 4, crosstalk = 0.15, gain = 1000, noise = 0.005, powerline_hz = 50, powerline_amp = 0.05, cmrr_db = 80, hp_cutoff = 10, offset = 0, bandwidth = c(20, 450), fs = 1000 )emgArray( channels = 4, crosstalk = 0.15, gain = 1000, noise = 0.005, powerline_hz = 50, powerline_amp = 0.05, cmrr_db = 80, hp_cutoff = 10, offset = 0, bandwidth = c(20, 450), fs = 1000 )
channels |
Number of electrode channels. |
crosstalk |
Per-neighbour crosstalk fraction (see |
gain |
Amplifier gain. |
noise |
Thermal-noise SD (referred to the electrode). |
powerline_hz, powerline_amp
|
Mains frequency (Hz) and its (common-mode) amplitude. |
cmrr_db |
Common-mode rejection ratio (dB) attenuating the mains interference. |
hp_cutoff |
DC-blocking high-pass corner (Hz). |
offset |
Electrode DC offset (a constant the DC-blocking high-pass removes). |
bandwidth |
Passband |
fs |
Sampling rate (Hz). |
an emg_array object.
emgArrayMeasure(), crosstalkMatrix(), emgElectrode()
emgArray(channels = 4, crosstalk = 0.2, cmrr_db = 90)emgArray(channels = 4, crosstalk = 0.2, cmrr_db = 90)
Turns a matrix of per-muscle activation envelopes into the multi-channel raw sEMG an electrode array records: each channel's activation-modulated Gaussian sEMG is band-limited, the channels are MIXED through the crosstalk matrix, and the amplifier adds CMRR-rejected mains interference, a DC-blocking high-pass, thermal noise and gain.
emgArrayMeasure(activations, sensor = emgArray(), seed = NULL)emgArrayMeasure(activations, sensor = emgArray(), seed = NULL)
activations |
An |
sensor |
An |
seed |
Optional RNG seed. |
an emg_array_measurement: raw (n_time x channels), clean
(uncontaminated, un-mixed, amplified), time, fs, the crosstalk matrix, sensor.
act <- matrix(pmax(0, sin(outer(seq(0, 6, len = 1200), 1:3))), ncol = 3) m <- emgArrayMeasure(act, emgArray(channels = 3, crosstalk = 0.2), seed = 1) dim(m$raw)act <- matrix(pmax(0, sin(outer(seq(0, 6, len = 1200), 1:3))), ncol = 3) m <- emgArrayMeasure(act, emgArray(channels = 3, crosstalk = 0.2), seed = 1) dim(m$raw)
Bundles the surface-EMG measurement chain: passband, powerline interference, baseline/motion wander, thermal noise, amplifier gain and clipping.
emgElectrode( gain = 1000, noise = 0.005, powerline_hz = 50, powerline_amp = 0.05, baseline_amp = 0.02, baseline_hz = 1, bandwidth = c(20, 450), fs = 1000, clip = Inf )emgElectrode( gain = 1000, noise = 0.005, powerline_hz = 50, powerline_amp = 0.05, baseline_amp = 0.02, baseline_hz = 1, bandwidth = c(20, 450), fs = 1000, clip = Inf )
gain |
Amplifier gain. |
noise |
Thermal (white) noise SD, referred to the electrode. |
powerline_hz |
Mains frequency (Hz); odd harmonics are added. |
powerline_amp |
Powerline interference amplitude. |
baseline_amp, baseline_hz
|
Baseline/motion-wander amplitude and frequency (Hz). |
bandwidth |
Passband |
fs |
Sampling rate (Hz). |
clip |
Saturation level (|output| clipped to this; |
an emg_electrode object.
emgMeasure(), emgEnvelope(), imuSensor()
emgElectrode(powerline_hz = 60, powerline_amp = 0.08)emgElectrode(powerline_hz = 60, powerline_amp = 0.08)
The standard sEMG envelope: full-wave rectification followed by a low-pass at
cutoff Hz (FFT-based), recovering the activation.
emgEnvelope(x, fs, cutoff = 5)emgEnvelope(x, fs, cutoff = 5)
x |
A raw or clean EMG signal. |
fs |
Sampling rate (Hz). |
cutoff |
Low-pass cut-off (Hz). |
the linear envelope (same length as x).
emgEnvelope(rnorm(1000) * abs(sin(seq(0, 6, length.out = 1000))), 1000)emgEnvelope(rnorm(1000) * abs(sin(seq(0, 6, length.out = 1000))), 1000)
Turns an activation envelope into a realistic raw sEMG: an activation-modulated Gaussian process band-passed to the electrode's passband (the clean EMG), then contaminated with powerline interference, baseline/motion wander and thermal noise and amplified.
emgMeasure(activation, sensor = emgElectrode(), seed = NULL)emgMeasure(activation, sensor = emgElectrode(), seed = NULL)
activation |
Non-negative activation envelope (0..1), length |
sensor |
An |
seed |
Optional RNG seed. |
an emg_measurement: raw (contaminated, amplified sEMG), clean
(uncontaminated band-limited EMG, amplified), time, activation, fs, sensor.
act <- pmax(0, sin(seq(0, 4 * pi, length.out = 2000))) emg <- emgMeasure(act, emgElectrode(fs = 1000), seed = 1) cor(emgEnvelope(emg$clean, 1000), act)act <- pmax(0, sin(seq(0, 4 * pi, length.out = 2000))) emg <- emgMeasure(act, emgElectrode(fs = 1000), seed = 1) cor(emgEnvelope(emg$clean, 1000), act)
A stochastic ensemble Kalman filter for x_k = f(x_{k-1}, k) + w,
y_k = h(x_k, k) + v. An ensemble of members carries the state distribution;
its sample cross- and innovation-covariances form the Kalman gain, and the
members are updated with perturbed observations. Scales to high-dimensional
states where a full covariance is impractical. Same interface as
unscentedKalmanFilter().
ensembleKalmanFilter(f, h, y, Q, R, x0, P0, n_ensemble = 100, seed = NULL)ensembleKalmanFilter(f, h, y, Q, R, x0, P0, n_ensemble = 100, seed = NULL)
f, h
|
Process and measurement models |
y |
Measurement matrix ( |
Q, R
|
Process- and measurement-noise covariances. |
x0, P0
|
Initial state mean and covariance. |
n_ensemble |
Ensemble size. |
seed |
Optional RNG seed. |
an enkf_result: state (n_steps x L ensemble means), cov (final
ensemble covariance), innovation.
Evensen G (1994) J Geophys Res 99:10143-10162; Burgers et al. (1998).
unscentedKalmanFilter(), extendedKalmanFilter(), particleFilter()
set.seed(1); y <- matrix(sin(seq(0, 6, length.out = 50)) + rnorm(50, 0, .1), 50, 1) fit <- ensembleKalmanFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1), n_ensemble = 80) fit$state[50]set.seed(1); y <- matrix(sin(seq(0, 6, length.out = 50)) + rnorm(50, 0, .1), 50, 1) fit <- ensembleKalmanFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1), n_ensemble = 80) fit$state[50]
After an intervention is delivered and a follow-up outcome is measured, check whether the observed change fell within the twin's predicted interval – the empirical test that closes the personalise -> predict -> act -> evaluate loop and calibrates trust in the twin.
evaluateIntervention( twin, chosen, observed_delta, readiness = NULL, outcome = "rom", outcome_sd = 0, n_mc = 200, duration = 6, dt = 0.005, theta0 = NULL, seed = NULL )evaluateIntervention( twin, chosen, observed_delta, readiness = NULL, outcome = "rom", outcome_sd = 0, n_mc = 200, duration = 6, dt = 0.005, theta0 = NULL, seed = NULL )
twin |
The personalised |
chosen |
The delivered |
observed_delta |
The measured outcome change at follow-up. |
readiness |
A |
outcome, n_mc, duration, dt, theta0, seed
|
As in |
outcome_sd |
Outcome observation/session variability (SD). The check uses a PREDICTION interval combining the parameter (epistemic) uncertainty with this (aleatoric) noise, so a realistic follow-up is covered – a parameter-only interval is over-narrow for a real measurement. |
an intervention_evaluation: predicted effect and prediction ci, the
observed change, and whether it was within_interval.
scoreInterventions(), clinicalDecision()
truth <- limbTwin(strength = 6) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01) rdy <- twinReadiness(fit, imu, dt = 0.01) evaluateIntervention(fit$twin, defaultInterventions()$strengthening, 1.5, rdy, outcome_sd = 0.1, n_mc = 50)truth <- limbTwin(strength = 6) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01) rdy <- twinReadiness(fit, imu, dt = 0.01) evaluateIntervention(fit$twin, defaultInterventions()$strengthening, 1.5, rdy, outcome_sd = 0.1, n_mc = 50)
A first-order extended Kalman filter for x_k = f(x_{k-1}, k) + w,
y_k = h(x_k, k) + v. The process and measurement Jacobians are obtained by
central finite differences unless analytic versions are supplied. Same
interface and return shape as unscentedKalmanFilter(); for a linear-Gaussian
model it reproduces the exact Kalman filter.
extendedKalmanFilter(f, h, y, Q, R, x0, P0, Fjac = NULL, Hjac = NULL)extendedKalmanFilter(f, h, y, Q, R, x0, P0, Fjac = NULL, Hjac = NULL)
f, h
|
Process and measurement models |
y |
Measurement matrix ( |
Q, R
|
Process- and measurement-noise covariances. |
x0, P0
|
Initial state mean and covariance. |
Fjac, Hjac
|
Optional analytic Jacobians |
an ekf_result: state (n_steps x L), cov (final), innovation.
Jazwinski AH (1970) Stochastic Processes and Filtering Theory.
unscentedKalmanFilter(), particleFilter(), ensembleKalmanFilter()
set.seed(1); y <- matrix(cumsum(rnorm(40, 0, .1)) + rnorm(40, 0, .1), 40, 1) fit <- extendedKalmanFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1)) fit$state[40]set.seed(1); y <- matrix(cumsum(rnorm(40, 0, .1)) + rnorm(40, 0, .1), 40, 1) fit <- extendedKalmanFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1)) fit$state[40]
Linear-interpolates the NA runs (occlusions) of a markerMeasure() result,
column by column, holding the endpoints flat where a gap reaches the boundary.
fillMarkerGaps(m)fillMarkerGaps(m)
m |
a |
the trajectory matrix with gaps interpolated.
pos <- matrix(cumsum(rnorm(200)), 200, 1) filled <- fillMarkerGaps(markerMeasure(pos, markerSensor(dropout_prob = 0.05), seed = 1))pos <- matrix(cumsum(rnorm(200)), 200, 1) filled <- fillMarkerGaps(markerMeasure(pos, markerSensor(dropout_prob = 0.05), seed = 1))
Samples n virtual subjects: physical parameters are drawn from the given
ranges, the IMU sensor is randomised (domain randomisation for robustness),
each twin is simulated and measured, and the sensor features are returned with
the generating parameters as labels. Use it to build ML training data where
the labels are true mechanistic quantities.
generateTrainingData( n = 100, strength = c(3, 10), damping = c(0.1, 0.8), stiffness = c(1.5, 4), duration = 5, dt = 0.02, randomize_sensor = TRUE, features = TRUE, seed = NULL )generateTrainingData( n = 100, strength = c(3, 10), damping = c(0.1, 0.8), stiffness = c(1.5, 4), duration = 5, dt = 0.02, randomize_sensor = TRUE, features = TRUE, seed = NULL )
n |
Number of virtual subjects. |
strength, damping, stiffness
|
Length-2 ranges (uniform) for each parameter. |
duration, dt
|
Simulation settings. |
randomize_sensor |
If |
features |
If |
seed |
Optional RNG seed. |
a twin_dataset: a data frame of features + label columns
(strength, damping, stiffness), or (when features = FALSE) a list
with labels and streams.
d <- generateTrainingData(20, duration = 4, dt = 0.02, seed = 1) head(d)d <- generateTrainingData(20, duration = 4, dt = 0.02, seed = 1) head(d)
Minimises an expensive loss over a bounded parameter box by Bayesian
optimisation: it fits a gpEmulator() surrogate to the evaluated
(parameters, loss) points and repeatedly evaluates the loss at the point of
maximum expected improvement. Far fewer loss evaluations than grid search or a
local optimiser – the tool for calibrating a slow simulator to data.
gpCalibrate( loss, ranges, n_init = 12, n_iter = 20, xi = 0.01, n_cand = 500, seed = NULL )gpCalibrate( loss, ranges, n_init = 12, n_iter = 20, xi = 0.01, n_cand = 500, seed = NULL )
loss |
|
ranges |
Named list of |
n_init |
Number of initial random design points. |
n_iter |
Number of Bayesian-optimisation iterations (loss evaluations after the initial design). |
xi |
Expected-improvement exploration margin. |
n_cand |
Candidate points scored by expected improvement each iteration. |
seed |
Optional RNG seed. |
a gp_calibration: the best par and value, the full design
X/y, the best-so-far trace, and the final gp surrogate.
Jones DR, Schonlau M, Welch WJ (1998) J Glob Optim 13:455-492 (efficient global optimisation / expected improvement).
gpEmulator(), abcSMC(), personalizeCardio()
# recover the minimum of a 2-D loss in a few evaluations loss <- function(p) (p[1] - 0.3)^2 + (p[2] - 0.7)^2 cal <- gpCalibrate(loss, list(x = c(0, 1), y = c(0, 1)), n_init = 8, n_iter = 15, seed = 1) round(cal$par, 2)# recover the minimum of a 2-D loss in a few evaluations loss <- function(p) (p[1] - 0.3)^2 + (p[2] - 0.7)^2 cal <- gpCalibrate(loss, list(x = c(0, 1), y = c(0, 1)), n_init = 8, n_iter = 15, seed = 1) round(cal$par, 2)
Fits a Gaussian-process regression surrogate with a squared-exponential kernel
to training inputs X and outputs y – typically a design of simulator runs,
so the emulator stands in for an expensive twin. Hyperparameters (length-scale,
signal and noise SD) are fitted by maximising the marginal likelihood.
predict.gpEmulator() returns the posterior mean and standard deviation at new
inputs; the SD is small near training points and grows away from them.
gpEmulator( X, y, optimize = TRUE, lengthscale = NULL, sigma_f = NULL, sigma_n = NULL )gpEmulator( X, y, optimize = TRUE, lengthscale = NULL, sigma_f = NULL, sigma_n = NULL )
X |
Training inputs: a matrix ( |
y |
Training outputs (length |
optimize |
If |
lengthscale, sigma_f, sigma_n
|
Initial (or fixed) kernel length-scale,
signal SD and noise SD; |
a gpEmulator object (with a predict.gpEmulator() method): the fitted
lengthscale, sigma_f, sigma_n and the training data.
Rasmussen CE, Williams CKI (2006) Gaussian Processes for Machine Learning.
predict.gpEmulator(), sobolIndices(), abcCalibration()
set.seed(1); x <- runif(20, 0, 6); y <- sin(x) + rnorm(20, 0, 0.05) gp <- gpEmulator(x, y) p <- predict(gp, seq(0, 6, length.out = 5)); round(p$mean, 2)set.seed(1); x <- runif(20, 0, 6); y <- sin(x) + rnorm(20, 0, 0.05) gp <- gpEmulator(x, y) p <- predict(gp, seq(0, 6, length.out = 5)); round(p$mean, 2)
Turns a simulated twin trajectory into a realistic IMU data stream by applying the forward sensor model and the sensor's imperfections.
imuMeasure(trajectory, sensor = imuSensor(), seed = NULL, noise = TRUE)imuMeasure(trajectory, sensor = imuSensor(), seed = NULL, noise = TRUE)
trajectory |
A |
sensor |
An |
seed |
Optional RNG seed. |
noise |
If |
a data frame time, gyro, ax, ay (angular velocity rad/s, body-frame
accelerations m/s^2).
imuSensor(), personalizeTwin()
sim <- simulateTwin(limbTwin(), duration = 2, dt = 0.01) head(imuMeasure(sim, imuSensor(gyro_noise = 0.04), seed = 1))sim <- simulateTwin(limbTwin(), duration = 2, dt = 0.01) head(imuMeasure(sim, imuSensor(gyro_noise = 0.04), seed = 1))
Bundles the imperfections of an inertial measurement unit – bias, white noise,
random-walk drift and scale error – so the same sensor can be reused, and
randomised across a population (see generateTrainingData()).
imuSensor( gyro_bias = 0, gyro_noise = 0.03, gyro_drift = 0, gyro_scale = 1, accel_bias = 0, accel_noise = 0.15, accel_scale = 1 )imuSensor( gyro_bias = 0, gyro_noise = 0.03, gyro_drift = 0, gyro_scale = 1, accel_bias = 0, accel_noise = 0.15, accel_scale = 1 )
gyro_bias, gyro_noise, gyro_drift, gyro_scale
|
Gyroscope constant bias (rad/s), white-noise SD (rad/s), random-walk drift SD per step (rad/s), and multiplicative scale error (1 = perfect). |
accel_bias, accel_noise, accel_scale
|
Accelerometer bias (m/s^2, length 1 or 2), white-noise SD (m/s^2) and scale error. |
an imu_sensor object.
imuSensor(gyro_bias = 0.02, gyro_noise = 0.05)imuSensor(gyro_bias = 0.02, gyro_noise = 0.05)
Applies parameter changes to a (personalised) twin – multiplicatively via
scale and/or by absolute values via set – re-simulates, and reports how
the movement outcomes change. This is the digital twin used as an experiment:
"if this subject's muscle were 30% stronger, how would the movement change?".
insilicoIntervention( twin, scale = list(), set = list(), duration = 6, dt = 0.005, theta0 = 0.1, omega0 = 0 )insilicoIntervention( twin, scale = list(), set = list(), duration = 6, dt = 0.005, theta0 = 0.1, omega0 = 0 )
twin |
A |
scale |
Named list of MULTIPLICATIVE changes, e.g. |
set |
Named list of ABSOLUTE parameter values to set. |
duration, dt, theta0, omega0
|
Simulation settings (as |
an insilico_intervention: before, after (trajectories),
outcomes (before/after/relative-change table) and the modified twin.
personalizeTwin(), simulateTwin()
tw <- limbTwin(strength = 5, damping = 0.6, stiffness = 3) iv <- insilicoIntervention(tw, scale = list(strength = 1.3, damping = 0.7)) iv$outcomestw <- limbTwin(strength = 5, damping = 0.6, stiffness = 3) iv <- insilicoIntervention(tw, scale = list(strength = 1.3, damping = 0.7)) iv$outcomes
Applies a fitness-raising intervention (e.g. a rehabilitation / conditioning programme) to a multi-modal twin and predicts the coherent change across modalities: movement range of motion up, resting/mean heart rate down.
insilicoInterventionMultimodal(mt, fitness_gain = 0.2, ...)insilicoInterventionMultimodal(mt, fitness_gain = 0.2, ...)
mt |
a |
fitness_gain |
Increase in latent fitness (0..1, clamped). |
... |
Passed to |
an mm_intervention: the before/after cross-modal summaries and
their change.
multimodalTwin(), populationIntervention()
insilicoInterventionMultimodal(multimodalTwin(fitness = 0.4), fitness_gain = 0.3, seed = 1)$changeinsilicoInterventionMultimodal(multimodalTwin(fitness = 0.4), fitness_gain = 0.3, seed = 1)$change
A candidate intervention for the clinical loop: a named set of parameter
changes (multiplicative scale and/or absolute set, as insilicoIntervention())
plus its clinical burden and risk, and optionally an evidence-based effect.
intervention( name, scale = list(), set = list(), burden = 0, risk = 0, evidence = NULL )intervention( name, scale = list(), set = list(), burden = 0, risk = 0, evidence = NULL )
name |
Label. |
scale, set
|
Parameter changes (see |
burden, risk
|
Clinical burden and risk, each on a 0..1 scale. |
evidence |
Optional |
a twin_intervention.
defaultInterventions(), scoreInterventions()
intervention("strengthening", scale = list(strength = 1.4), burden = 0.4, risk = 0.1)intervention("strengthening", scale = list(strength = 1.4), burden = 0.4, risk = 0.1)
Simulates a single cortical column with the Jansen-Rit neural-mass model: three
interacting neural populations (pyramidal, excitatory and inhibitory
interneurons) whose average post-synaptic potentials evolve as second-order
("alpha-function") systems coupled through a sigmoidal rate-to-potential
nonlinearity, driven by a background firing-rate input p(t). With the default
parameters the column produces a self-organised ~10 Hz alpha rhythm; the
scalp EEG is taken as the pyramidal post-synaptic potential y1 - y2.
jansenRit( duration = 8, sfeeg = 256, A = 3.25, B = 22, a = 100, b = 50, C = 135, p_mean = 220, p_sd = 30, oversample = 4, aperiodic = 0, aperiodic_exp = 1.4, seed = NULL )jansenRit( duration = 8, sfeeg = 256, A = 3.25, B = 22, a = 100, b = 50, C = 135, p_mean = 220, p_sd = 30, oversample = 4, aperiodic = 0, aperiodic_exp = 1.4, seed = NULL )
duration |
Recording length (s). |
sfeeg |
Output sampling rate (Hz). |
A, B
|
Excitatory and inhibitory average synaptic gains (mV). |
a, b
|
Excitatory and inhibitory rate constants (inverse time constants, 1/s). |
C |
Global connectivity constant; the four intra-column connectivities are
|
p_mean, p_sd
|
Mean and SD of the background firing-rate input |
oversample |
Internal integration steps per output sample (accuracy
margin for the fast |
aperiodic |
Amplitude of an added 1/f background, as a multiple of the
oscillation's SD (0 = the pure neural-mass rhythm; a real scalp channel needs
|
aperiodic_exp |
Spectral slope of the 1/f background (see |
seed |
Optional RNG seed. |
The excitation/inhibition balance and the drive select the regime: reducing the
mean input p_mean below the oscillation threshold collapses the alpha limit
cycle to a low-amplitude resting fixed point, while changing the gains A/B
or the connectivity C shifts the rhythm.
a jansenRit object: time, eeg (the y1 - y2 scalp potential, mV),
dominant_hz (the dominant rhythm above 3 Hz, robust to the 1/f floor),
band (its conventional EEG band), sfeeg and params.
Jansen BH, Rit VG (1995) Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern 73:357-366. Grimbert F, Faugeras O (2006) Neural Comput 18:3052-3068.
ecgsyn(), kuramoto(), cpgMatsuoka()
eeg <- jansenRit(duration = 6, seed = 1) eeg$dominant_hz # ~10 Hz alpha rhythm eeg$band # "alpha"eeg <- jansenRit(duration = 6, seed = 1) eeg$dominant_hz # ~10 Hz alpha rhythm eeg$band # "alpha"
Simulates n phase oscillators with natural frequencies omega and all-to-all
coupling K (theta_i' = omega_i + K/n sum_j sin(theta_j - theta_i)). The
order parameter measures collective synchronisation: near 0 (incoherent) below
the critical coupling, near 1 (synchronised) above it.
kuramoto(n = 20, K = 1, omega = NULL, duration = 20, dt = 0.01, seed = NULL)kuramoto(n = 20, K = 1, omega = NULL, duration = 20, dt = 0.01, seed = NULL)
n |
Number of oscillators. |
K |
Coupling strength. |
omega |
Natural frequencies (length |
duration, dt
|
Simulation length and step (s). |
seed |
Optional RNG seed. |
a kuramoto object: time, phase (nt x n), order (order
parameter per step), signals (cos(phase)), mean_order.
Kuramoto Y (1975); Strogatz SH (2000) Physica D 143:1-20.
cpgMatsuoka(), ecgsyn(), jansenRit()
kuramoto(n = 15, K = 3, duration = 20, seed = 1)$mean_order # synchronised kuramoto(n = 15, K = 0, duration = 20, seed = 1)$mean_order # incoherentkuramoto(n = 15, K = 3, duration = 20, seed = 1)$mean_order # synchronised kuramoto(n = 15, K = 0, duration = 20, seed = 1)$mean_order # incoherent
Constructs a mechanistic single-joint limb model whose parameters map to
physical/clinical quantities (segment inertia, active muscle strength, passive
damping, joint stiffness). It is the object simulated, personalised and
intervened on by the rest of PhysioTwin.
limbTwin( inertia = 0.15, strength = 5, damping = 0.3, stiffness = 2.5, freq = 1, rest = 0, mass = 2.5, com = 0.2, g = 9.81, sensor_r = 0.25 )limbTwin( inertia = 0.15, strength = 5, damping = 0.3, stiffness = 2.5, freq = 1, rest = 0, mass = 2.5, com = 0.2, g = 9.81, sensor_r = 0.25 )
inertia |
Segment moment of inertia about the joint (kg m^2). |
strength |
Active-drive amplitude (N m) – the muscle strength. |
damping |
Passive viscous damping (N m s / rad). |
stiffness |
Joint stiffness toward |
freq |
Rhythm of the active drive (Hz). |
rest |
Neutral joint angle (rad). |
mass, com, g
|
Segment mass (kg), centre-of-mass distance from the joint
(m) and gravity (m/s^2); set |
sensor_r |
Distance of the IMU from the joint along the segment (m), used by the sensor layer. |
a movement_twin object (a named parameter list with class).
simulateTwin(), personalizeTwin(), insilicoIntervention()
tw <- limbTwin(strength = 6, damping = 0.4, stiffness = 3) sim <- simulateTwin(tw, duration = 4, dt = 0.005) range(sim$theta)tw <- limbTwin(strength = 6, damping = 0.4, stiffness = 3) sim <- simulateTwin(tw, duration = 4, dt = 0.005) range(sim$theta)
Applies the forward marker sensor model to a clean trajectory: adds Gaussian
noise, carves contiguous occlusion gaps (visible/occluded Markov process),
optionally quantises and injects ghost samples. Occluded frames are returned as
NA rows.
markerMeasure(position, sensor = markerSensor(), seed = NULL)markerMeasure(position, sensor = markerSensor(), seed = NULL)
position |
Clean marker trajectory: an |
sensor |
A |
seed |
Optional RNG seed. |
a marker_measurement: position (measured, NA rows where occluded),
visible (logical per frame), gap_fraction, and the sensor.
markerSensor(), fillMarkerGaps()
pos <- cbind(sin(seq(0, 6, length.out = 300)), cos(seq(0, 6, length.out = 300))) m <- markerMeasure(pos, markerSensor(dropout_prob = 0.02, dropout_len = 12), seed = 1) m$gap_fractionpos <- cbind(sin(seq(0, 6, length.out = 300)), cos(seq(0, 6, length.out = 300))) m <- markerMeasure(pos, markerSensor(dropout_prob = 0.02, dropout_len = 12), seed = 1) m$gap_fraction
Bundles the imperfections of an optical marker: positional noise, occlusion dropout (contiguous gaps), quantisation and ghost samples.
markerSensor( noise = 0.002, dropout_prob = 0.02, dropout_len = 10, ghost_prob = 0, quantization = 0 )markerSensor( noise = 0.002, dropout_prob = 0.02, dropout_len = 10, ghost_prob = 0, quantization = 0 )
noise |
Gaussian position-noise SD (same units as the trajectory, e.g. m). |
dropout_prob |
Per-frame probability of STARTING an occlusion while visible. |
dropout_len |
Mean occlusion length (frames); gaps are geometric. |
ghost_prob |
Per-frame probability of a spurious (ghost) sample. |
quantization |
Position resolution (round to this grid; 0 = off). |
a marker_sensor object.
markerMeasure(), fillMarkerGaps(), imuSensor()
markerSensor(noise = 0.001, dropout_prob = 0.03, dropout_len = 8)markerSensor(noise = 0.001, dropout_prob = 0.03, dropout_len = 8)
Samples from a log-posterior by random-walk Metropolis-Hastings. The Gaussian proposal covariance adapts to the empirical covariance of the chain during the burn-in (Haario et al. 2001), then is frozen, so correlated targets are sampled efficiently while the stationary distribution is preserved.
metropolis( logpost, init, n_iter = 6000, proposal_sd = 0.5, burn_in = 0.4, seed = NULL )metropolis( logpost, init, n_iter = 6000, proposal_sd = 0.5, burn_in = 0.4, seed = NULL )
logpost |
|
init |
Initial parameter vector (its names label the output). |
n_iter |
Number of iterations. |
proposal_sd |
Initial per-parameter proposal SD. |
burn_in |
Fraction discarded (and over which the proposal adapts). |
seed |
Optional RNG seed. |
an mcmc_result: posterior samples, acceptance rate, mean, sd,
and 2.5/50/97.5% quantiles.
Metropolis et al. (1953); Haario H, Saksman E, Tamminen J (2001) Bernoulli 7:223-242 (adaptive Metropolis).
abcSMC(), abcCalibration(), gpCalibrate()
# posterior of a normal mean: 20 obs, N(0,10) prior set.seed(1); y <- rnorm(20, 2, 1) lp <- function(m) sum(dnorm(y, m, 1, log = TRUE)) + dnorm(m, 0, 10, log = TRUE) fit <- metropolis(lp, init = c(mu = 0), n_iter = 4000, seed = 1) fit$mean# posterior of a normal mean: 20 obs, N(0,10) prior set.seed(1); y <- rnorm(20, 2, 1) lp <- function(m) sum(dnorm(y, m, 1, log = TRUE)) + dnorm(m, 0, 10, log = TRUE) fit <- metropolis(lp, init = c(mu = 0), n_iter = 4000, seed = 1) fit$mean
Couples a latent fitness and effort to both a limbTwin() movement model
(higher fitness -> more strength, less damping) and cardiorespiratory set-points
(higher fitness -> lower resting heart rate; higher effort -> higher heart rate
and breathing rate).
multimodalTwin( fitness = 0.5, effort = 0.5, base_strength = 5, base_damping = 0.5 )multimodalTwin( fitness = 0.5, effort = 0.5, base_strength = 5, base_damping = 0.5 )
fitness |
Latent fitness (0..1). |
effort |
Latent effort/activity level (0..1). |
base_strength, base_damping
|
Reference movement parameters at |
a multimodal_twin: the derived movement twin, resting_hr,
resp_rate and the latent fitness/effort.
simulateMultimodal(), insilicoInterventionMultimodal(), rsaTachogram()
mt <- multimodalTwin(fitness = 0.7) mt$resting_hrmt <- multimodalTwin(fitness = 0.7) mt$resting_hr
Chooses the experimental condition that will best identify a movement twin's parameters. For each candidate value of a design variable (e.g. the driving frequency), it forms the Fisher information matrix of the noise-weighted IMU observable about the parameters and scores the design by an optimality criterion – D-optimality (the log-determinant of the information, maximising overall information) or A-optimality (minimising the total Cramer-Rao variance). The most informative condition is the experiment worth running to pin the parameters down.
optimalDesign( twin, estimate = "strength", designs, design_var = "freq", sensor = imuSensor(), duration = 6, dt = 0.01, theta0 = 0.1, criterion = c("D", "A") )optimalDesign( twin, estimate = "strength", designs, design_var = "freq", sensor = imuSensor(), duration = 6, dt = 0.01, theta0 = 0.1, criterion = c("D", "A") )
twin |
A |
estimate |
Parameter(s) whose identifiability is optimised. |
designs |
Numeric vector of candidate values for the design variable. |
design_var |
The twin field the design varies (e.g. |
sensor |
An |
duration, dt, theta0
|
Simulation settings for the observable. |
criterion |
|
an optimal_design: a table (design, logdet_fim, score, and a
crlb_<param> Cramer-Rao SD per estimated parameter), the best design, and
the estimate/design_var/criterion.
Chaloner K, Verdinelli I (1995) Bayesian experimental design: a review. Statist Sci 10:273-304; Fedorov (1972) D-optimality.
personalizeTwin(), validateTwin(), profileLikelihood()
od <- optimalDesign(limbTwin(), estimate = "strength", designs = seq(0.4, 2.4, by = 0.2), design_var = "freq", duration = 5, dt = 0.02) od$best # the most informative driving frequencyod <- optimalDesign(limbTwin(), estimate = "strength", designs = seq(0.4, 2.4, by = 0.2), design_var = "freq", duration = 5, dt = 0.02) od$best # the most informative driving frequency
A sequential-importance-resampling (bootstrap) particle filter for
x_k = f(x_{k-1}, k) + w, y_k = h(x_k, k) + v. Particles are propagated
through f with process noise Q, weighted by the Gaussian measurement
likelihood, and resampled (systematic) when the effective sample size falls
below resample_thresh of the population. Handles strong nonlinearity and
multi-modal posteriors that a Kalman filter cannot. Same interface as
unscentedKalmanFilter().
particleFilter( f, h, y, Q, R, x0, P0, n_particles = 1000, resample_thresh = 0.5, seed = NULL )particleFilter( f, h, y, Q, R, x0, P0, n_particles = 1000, resample_thresh = 0.5, seed = NULL )
f, h
|
Process and measurement models |
y |
Measurement matrix ( |
Q, R
|
Process- and measurement-noise covariances. |
x0, P0
|
Initial state mean and covariance. |
n_particles |
Number of particles. |
resample_thresh |
Resample when effective sample size |
seed |
Optional RNG seed. |
a pf_result: state (n_steps x L posterior means), cov (final
weighted covariance), innovation, and ess (effective sample size per step).
Gordon NJ, Salmond DJ, Smith AFM (1993) IEE Proc F 140:107-113.
unscentedKalmanFilter(), extendedKalmanFilter(), ensembleKalmanFilter()
set.seed(1); y <- matrix(sin(seq(0, 6, length.out = 50)) + rnorm(50, 0, .1), 50, 1) fit <- particleFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1), n_particles = 500) fit$state[50]set.seed(1); y <- matrix(sin(seq(0, 6, length.out = 50)) + rnorm(50, 0, .1), 50, 1) fit <- particleFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1), n_particles = 500) fit$state[50]
Samples the joint posterior of a state-space model's latent trajectory and parameters by particle Gibbs: each iteration draws a new state trajectory with a conditional SMC (ancestor-sampling) sweep given the current parameters, then moves the parameters with a Metropolis step using the tractable complete-data likelihood of that trajectory. Ancestor sampling gives good mixing even for long series where plain particle Gibbs degenerates.
particleGibbs( build_model, y, prior_logd, init, n_particles = 100, n_iter = 2000, proposal_sd = 0.1, burn_in = 0.3, seed = NULL )particleGibbs( build_model, y, prior_logd, init, n_particles = 100, n_iter = 2000, proposal_sd = 0.1, burn_in = 0.3, seed = NULL )
build_model |
|
y |
Observation matrix ( |
prior_logd |
|
init |
Initial parameter vector (its names label the output). |
n_particles |
Particles in the conditional SMC. |
n_iter |
Iterations. |
proposal_sd |
Initial per-parameter proposal SD. |
burn_in |
Fraction discarded (and over which the proposal adapts). |
seed |
Optional RNG seed. |
a pgibbs_result: parameter posterior samples, acceptance, mean,
sd, 2.5/50/97.5% quantiles, and the posterior-mean state trajectory.
Andrieu C, Doucet A, Holenstein R (2010) J R Stat Soc B 72:269-342; Lindsten F, Jordan MI, Schon TB (2014) J Mach Learn Res 15:2145-2184.
particleMCMC(), metropolis(), particleFilter()
set.seed(1); a <- 0.8; x <- numeric(50); for (k in 2:50) x[k] <- a * x[k-1] + rnorm(1, 0, 0.3) y <- matrix(x + rnorm(50, 0, 0.3), 50, 1) bm <- function(th) list(f = function(x, k) th[1] * x, h = function(x, k) x, Q = matrix(0.09), R = matrix(0.09), x0 = 0, P0 = matrix(1)) fit <- particleGibbs(bm, y, function(th) dunif(th[1], -1, 1, log = TRUE), init = c(a = 0.4), n_particles = 60, n_iter = 600, seed = 1) fit$meanset.seed(1); a <- 0.8; x <- numeric(50); for (k in 2:50) x[k] <- a * x[k-1] + rnorm(1, 0, 0.3) y <- matrix(x + rnorm(50, 0, 0.3), 50, 1) bm <- function(th) list(f = function(x, k) th[1] * x, h = function(x, k) x, Q = matrix(0.09), R = matrix(0.09), x0 = 0, P0 = matrix(1)) fit <- particleGibbs(bm, y, function(th) dunif(th[1], -1, 1, log = TRUE), init = c(a = 0.4), n_particles = 60, n_iter = 600, seed = 1) fit$mean
Samples the parameter posterior of a state-space model from a WHOLE observation
series, when the likelihood is intractable. Each proposed parameter is scored by
a bootstrap particle filter's unbiased estimate of the marginal likelihood
p(y_{1:T} | theta), which enters the Metropolis-Hastings acceptance ratio;
the resulting chain targets the exact parameter posterior (Andrieu-Doucet-
Holenstein). The proposal covariance adapts during burn-in.
particleMCMC( build_model, y, prior_logd, init, n_particles = 200, n_iter = 3000, proposal_sd = 0.1, burn_in = 0.3, seed = NULL )particleMCMC( build_model, y, prior_logd, init, n_particles = 200, n_iter = 3000, proposal_sd = 0.1, burn_in = 0.3, seed = NULL )
build_model |
|
y |
Observation matrix ( |
prior_logd |
|
init |
Initial parameter vector (its names label the output). |
n_particles |
Particles in the likelihood-estimating filter. |
n_iter |
MCMC iterations. |
proposal_sd |
Initial per-parameter proposal SD. |
burn_in |
Fraction discarded (and over which the proposal adapts). |
seed |
Optional RNG seed. |
a pmcmc_result: posterior samples, acceptance, mean, sd and
2.5/50/97.5% quantiles.
Andrieu C, Doucet A, Holenstein R (2010) J R Stat Soc B 72:269-342.
particleFilter(), metropolis(), personalizeCardio()
# recover the AR coefficient of a linear-Gaussian state-space model from the series set.seed(1); a <- 0.7; x <- numeric(60); for (k in 2:60) x[k] <- a * x[k-1] + rnorm(1, 0, 0.3) y <- matrix(x + rnorm(60, 0, 0.3), 60, 1) bm <- function(th) list(f = function(x, k) th[1] * x, h = function(x, k) x, Q = matrix(0.09), R = matrix(0.09), x0 = 0, P0 = matrix(1)) fit <- particleMCMC(bm, y, function(th) dunif(th[1], -1, 1, log = TRUE), init = c(a = 0.4), n_particles = 100, n_iter = 800, seed = 1) fit$mean# recover the AR coefficient of a linear-Gaussian state-space model from the series set.seed(1); a <- 0.7; x <- numeric(60); for (k in 2:60) x[k] <- a * x[k-1] + rnorm(1, 0, 0.3) y <- matrix(x + rnorm(60, 0, 0.3), 60, 1) bm <- function(th) list(f = function(x, k) th[1] * x, h = function(x, k) x, Q = matrix(0.09), R = matrix(0.09), x0 = 0, P0 = matrix(1)) fit <- particleMCMC(bm, y, function(th) dunif(th[1], -1, 1, log = TRUE), init = c(a = 0.4), n_particles = 100, n_iter = 800, seed = 1) fit$mean
Fits the closed-loop cardioRespiratory() model's baroreflex gains and
respiratory modulation to a subject's measured variability summaries (baroreflex
sensitivity, blood-pressure Mayer power, heart-rate respiratory power) by
GP-accelerated calibration (gpCalibrate()). Because the summaries have very
different scales, the discrepancy is relative. The returned twin, driven with
the fitted parameters, reproduces the subject's HRV/BPV.
personalizeCardio( target, ranges = list(g_vagal = c(0.004, 0.02), g_resist = c(0.003, 0.02), rsa_amp = c(0.01, 0.06)), fixed = list(), duration = 300, n_init = 12, n_iter = 25, n_rep = 3, seed = 1 )personalizeCardio( target, ranges = list(g_vagal = c(0.004, 0.02), g_resist = c(0.003, 0.02), rsa_amp = c(0.01, 0.06)), fixed = list(), duration = 300, n_init = 12, n_iter = 25, n_rep = 3, seed = 1 )
target |
Named numeric vector of the observed summaries to match; any of
|
ranges |
Named list of |
fixed |
Named list of other |
duration |
Simulation length per evaluation (s). |
n_init, n_iter
|
Bayesian-optimisation initial design and iterations. |
n_rep |
Simulations averaged per evaluation to stabilise the (stochastic) variability summaries – the blood-pressure Mayer power is noisy, so averaging makes the inverse problem well-posed. |
seed |
RNG seed (fixes the stochastic loss so the calibration is deterministic). |
a cardio_personalization: the fitted estimates, the personalised
twin, the target and achieved summaries, the relative error rel_error,
and the calibration object.
cardioRespiratory(), cardioSummary(), gpCalibrate(), personalizeTwin()
subject <- cardioSummary(cardioRespiratory(duration = 240, g_vagal = 0.013, g_resist = 0.008, rsa_amp = 0.045, seed = 3)) fit <- personalizeCardio(subject, n_init = 8, n_iter = 15, seed = 1) fit$achieved # reproduces the subject's brs / bp_lf / hr_hfsubject <- cardioSummary(cardioRespiratory(duration = 240, g_vagal = 0.013, g_resist = 0.008, rsa_amp = 0.045, seed = 3)) fit <- personalizeCardio(subject, n_init = 8, n_iter = 15, seed = 1) fit$achieved # reproduces the subject's brs / bp_lf / hr_hf
Fits the cardioStateSpace() parameters jointly to a whole RR-interval record.
Because the closed-loop model is linear-Gaussian its likelihood is exact (a
Kalman filter), so the posterior is sampled by metropolis() – the
whole-waveform counterpart to the summary-based personalizeCardio(). The Mayer
gain/damping and the respiratory amplitude are identifiable here because the
record separates the low- and high-frequency rhythms, resolving the ridge that
confounds the summary fit. (For an intractable, nonlinear model use
particleMCMC() / particleGibbs() with cardioStateSpace()'s f/h.)
personalizeCardioWaveform( rr, resp_freq = 0.25, estimate = NULL, lf_model = c("oscillator", "ar1"), fixed = list(omega_m = 0.53, rr0 = NULL, q_mayer = 1, r_obs = 9), n_iter = 4000, seed = NULL )personalizeCardioWaveform( rr, resp_freq = 0.25, estimate = NULL, lf_model = c("oscillator", "ar1"), fixed = list(omega_m = 0.53, rr0 = NULL, q_mayer = 1, r_obs = 9), n_iter = 4000, seed = NULL )
rr |
RR-interval series (ms). |
resp_freq |
Respiratory frequency (Hz) for the RSA regressor. |
estimate |
Parameters to fit (any of |
lf_model |
Low-frequency model (see |
fixed |
Named list of held-fixed |
n_iter |
Metropolis iterations. |
seed |
Optional RNG seed. |
a cardio_waveform_fit: the posterior (an mcmc_result), the fitted
estimates, and the estimated names.
cardioStateSpace(), metropolis(), particleMCMC(), personalizeCardio()
truth <- cardioStateSpace(g = 8, rho = 0.9, A_rsa = 20, resp = sin(2 * pi * 0.25 * cumsum(rep(0.85, 250)))) set.seed(1); x <- c(0, 0); rr <- numeric(250) for (k in seq_len(250)) { x <- truth$f(x, k) + c(rnorm(1), 0); rr[k] <- truth$h(x, k) + rnorm(1, 0, 3) } fit <- personalizeCardioWaveform(rr, resp_freq = 0.25, estimate = c("g", "A_rsa"), n_iter = 3000, seed = 1) fit$estimatestruth <- cardioStateSpace(g = 8, rho = 0.9, A_rsa = 20, resp = sin(2 * pi * 0.25 * cumsum(rep(0.85, 250)))) set.seed(1); x <- c(0, 0); rr <- numeric(250) for (k in seq_len(250)) { x <- truth$f(x, k) + c(rnorm(1), 0); rr[k] <- truth$h(x, k) + rnorm(1, 0, 3) } fit <- personalizeCardioWaveform(rr, resp_freq = 0.25, estimate = c("g", "A_rsa"), n_iter = 3000, seed = 1) fit$estimates
Estimates a subject's twin parameters by assimilating their IMU stream with an
unscented Kalman filter over an augmented state [theta, omega, log(params)].
Parameters are estimated in log-space (keeping them positive) and the fitted
values are read from the converged filter.
personalizeTwin( imu, prior, estimate = "strength", method = c("optim", "ukf"), dt = 0.005, param_rw = 0.02, gyro_var = 0.03^2, accel_var = 0.15^2, burn_in = 0.5, n_starts = 25L )personalizeTwin( imu, prior, estimate = "strength", method = c("optim", "ukf"), dt = 0.005, param_rw = 0.02, gyro_var = 0.03^2, accel_var = 0.15^2, burn_in = 0.5, n_starts = 25L )
imu |
A data frame with |
prior |
A |
estimate |
Character vector of parameters to personalise (any of
|
method |
|
dt |
Sample period (s). |
param_rw |
Random-walk SD (log-space) for the estimated parameters (their assumed drift; default 0.02). |
gyro_var, accel_var
|
Measurement-noise variances for the gyro and accel channels. |
burn_in |
Fraction of the record discarded before averaging the converged estimate (default 0.5). |
n_starts |
Number of coarse-grid starting points scanned before the
|
a personalized_twin: twin (the fitted twin), estimates (named
fitted parameters), the co-estimated gyro_bias, and – for method = "ukf" – the per-step trajectory and the ukf result, or – for method = "optim" – the fitted theta0 and final objective.
unscentedKalmanFilter(), extendedKalmanFilter(), particleFilter(),
abcCalibration(), insilicoIntervention()
truth <- limbTwin(strength = 7, damping = 0.5) sim <- simulateTwin(truth, duration = 6, dt = 0.01) imu <- imuMeasure(sim, imuSensor(gyro_noise = 0.03, accel_noise = 0.1), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4, damping = 0.2), estimate = c("strength", "damping"), dt = 0.01) fit$estimatestruth <- limbTwin(strength = 7, damping = 0.5) sim <- simulateTwin(truth, duration = 6, dt = 0.01) imu <- imuMeasure(sim, imuSensor(gyro_noise = 0.03, accel_noise = 0.1), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4, damping = 0.2), estimate = c("strength", "damping"), dt = 0.01) fit$estimates
Samples a cohort of movement twins from a parameter distribution, applies the
same intervention to each (via insilicoIntervention()), and summarises the
response as a paired in-silico trial: the mean effect and its paired effect size
(Cohen's d_z), the responder rate against a minimal clinically important
difference, the statistical power at this sample size and the sample size needed
for 80% power.
populationIntervention( n = 200, params = list(strength = c(5, 1), damping = c(0.7, 0.2), stiffness = c(2.5, 0.6)), intervention = list(scale = list(strength = 1.4)), outcome = "rom", mcid = NULL, duration = 6, dt = 0.005, theta0 = 0.1, omega0 = 0, seed = NULL )populationIntervention( n = 200, params = list(strength = c(5, 1), damping = c(0.7, 0.2), stiffness = c(2.5, 0.6)), intervention = list(scale = list(strength = 1.4)), outcome = "rom", mcid = NULL, duration = 6, dt = 0.005, theta0 = 0.1, omega0 = 0, seed = NULL )
n |
Number of subjects (twins) in the cohort. |
params |
Named list of |
intervention |
List with |
outcome |
Outcome to track: one of |
mcid |
Minimal clinically important difference (responder threshold); defaults to 0.1 x the baseline outcome SD. |
duration, dt, theta0, omega0
|
Simulation settings (as |
seed |
Optional RNG seed. |
a population_intervention: subjects (per-subject params, pre,
post, delta), effect (mean delta, SD, Cohen's d_z), responder_rate,
mcid, power and n_for_80.
insilicoIntervention(), multimodalTwin(), generateTrainingData()
pop <- populationIntervention(n = 100, intervention = list(scale = list(strength = 1.4)), seed = 1) pop$effectpop <- populationIntervention(n = 100, intervention = list(scale = list(strength = 1.4)), seed = 1) pop$effect
Predict from a Gaussian-process emulator
## S3 method for class 'gpEmulator' predict(object, newdata, ...)## S3 method for class 'gpEmulator' predict(object, newdata, ...)
object |
a |
newdata |
New inputs: a matrix ( |
... |
Unused. |
a list with mean (posterior mean) and sd (posterior standard
deviation) at each new input.
Assesses whether a parameter is practically identifiable from a recording:
it is fixed across a grid of values and the remaining nuisance parameters
(plus a start angle and gyro bias) are re-optimised, tracing the fit cost. A
clear minimum bounded by the chi-square threshold gives a confidence interval;
a flat profile means the parameter cannot be pinned down.
profileLikelihood( imu, prior, param, over = c(0.4, 2), nuisance = character(0), n_grid = 15L, dt = 0.02, gyro_var = 0.03^2, accel_var = 0.12^2 )profileLikelihood( imu, prior, param, over = c(0.4, 2), nuisance = character(0), n_grid = 15L, dt = 0.02, gyro_var = 0.03^2, accel_var = 0.12^2 )
imu |
An IMU recording ( |
prior |
A |
param |
The parameter to profile. |
over |
Length-2 multiplicative range around |
nuisance |
Other parameters re-optimised at each grid point. |
n_grid |
Grid size (default 15). |
dt, gyro_var, accel_var
|
Sample period and measurement variances. |
a profile_likelihood: the profile (grid value, cost), the cost
minimum, the 95% ci (or NA if open), and identifiable (logical).
Raue A, et al. (2009) Bioinformatics 25:1923-1929.
validateTwin(), sobolIndices()
tw <- limbTwin(strength = 6, damping = 0.5, stiffness = 6) imu <- imuMeasure(simulateTwin(tw, 8, 0.02, theta0 = 0.2), imuSensor(gyro_noise = 0.03, accel_noise = 0.12), seed = 1) profileLikelihood(imu, tw, "strength", nuisance = "damping", n_grid = 9)$citw <- limbTwin(strength = 6, damping = 0.5, stiffness = 6) imu <- imuMeasure(simulateTwin(tw, 8, 0.02, theta0 = 0.2), imuSensor(gyro_noise = 0.03, accel_noise = 0.12), seed = 1) profileLikelihood(imu, tw, "strength", nuisance = "damping", n_grid = 9)$ci
The theoretical quantisation noise (least-significant-bit squared over twelve)
and the empirical error of quantising signal to bits over range.
quantizationNoise(signal, bits, range = NULL)quantizationNoise(signal, bits, range = NULL)
signal |
Numeric signal. |
bits |
Bit depth. |
range |
|
a named numeric vector: lsb, theoretical_sd (= lsb/sqrt(12)),
empirical_sd, theoretical_power (= lsb^2/12).
quantizationNoise(sin(seq(0, 20, length.out = 4000)), bits = 12)quantizationNoise(sin(seq(0, 20, length.out = 4000)), bits = 12)
Combines the twin's mechanistic effect estimate with an external (e.g. literature) estimate by inverse-variance (precision) weighting – the Bayesian conjugate-normal update – and flags whether the two AGREE. Evidence complements the mechanistic prediction; it is never overridden.
reconcileEvidence(mech_est, mech_se, ev_est, ev_se)reconcileEvidence(mech_est, mech_se, ev_est, ev_se)
mech_est, mech_se
|
Mechanistic effect estimate and its SE. |
ev_est, ev_se
|
Evidence-based effect estimate and its SE. |
an evidence_reconciliation: combined estimate/se, the discrepancy
z, and agreement (agree/disagree, at |z| < 2).
scoreInterventions(), abcCalibration()
reconcileEvidence(mech_est = 0.30, mech_se = 0.05, ev_est = 0.26, ev_se = 0.04)reconcileEvidence(mech_est = 0.30, mech_se = 0.05, ev_est = 0.26, ev_se = 0.04)
Connects the digital-twin decision loop to the PhysioRehab clinical-reasoning
framework. A real single-case series – outcome measurements in a baseline
phase then an intervention phase – is assembled into a rehab_episode and
analysed with PhysioRehab's single-case experimental design engine: the NAP and
Tau effect sizes, whether the change exceeds the mcid, and the ICF category of
the measure. If the twin's predicted_delta is supplied, it is reconciled
with the clinically-observed change (do they agree in sign, and on the MCID
verdict?), closing the predict -> deliver -> evaluate loop against real evidence.
rehabEvaluate( baseline, intervention, predicted_delta = NULL, measure = "gait_speed", mcid = 0.16, direction = "increase", condition = "in-silico twin", patient_id = "twin", start_date = "2026-01-01" )rehabEvaluate( baseline, intervention, predicted_delta = NULL, measure = "gait_speed", mcid = 0.16, direction = "increase", condition = "in-silico twin", patient_id = "twin", start_date = "2026-01-01" )
baseline |
Numeric outcome measurements in the baseline phase. |
intervention |
Numeric outcome measurements in the intervention phase. |
predicted_delta |
Optional: the twin's predicted change in the outcome
(e.g. from |
measure |
Outcome name (a PhysioRehab measure, e.g. |
mcid |
Minimal clinically important difference (threshold, outcome units). |
direction |
|
condition, patient_id
|
Episode metadata. |
start_date |
First session date ( |
a rehab_evaluation: the built episode, the PhysioRehab sced result,
the icf code(s), and – when predicted_delta is given – a reconciliation
(predicted, observed, sign_agree, mcid_agree, error).
evaluateIntervention(), clinicalDecision()
if (requireNamespace("PhysioRehab", quietly = TRUE)) { set.seed(1) base <- rnorm(5, 0.80, 0.03); post <- rnorm(6, 1.00, 0.03) rehabEvaluate(base, post, predicted_delta = 0.18, measure = "gait_speed", mcid = 0.16) }if (requireNamespace("PhysioRehab", quietly = TRUE)) { set.seed(1) base <- rnorm(5, 0.80, 0.03); post <- rnorm(6, 1.00, 0.03) rehabEvaluate(base, post, predicted_delta = 0.18, measure = "gait_speed", mcid = 0.16) }
Simulates a breathing (chest-volume) waveform: within each breath the volume
rises smoothly on inspiration and falls on expiration, with the inspiration
fraction set by the inspiration:expiration (ie_ratio) ratio (physiological
breathing is expiration-dominant, ~1:2).
respiration(duration = 60, rate = 15, ie_ratio = 0.5, amplitude = 1, sfp = 100)respiration(duration = 60, rate = 15, ie_ratio = 0.5, amplitude = 1, sfp = 100)
duration |
Recording length (s). |
rate |
Breathing rate (breaths/min). |
ie_ratio |
Inspiration-to-expiration duration ratio (0.5 = 1:2). |
amplitude |
Peak-to-trough volume amplitude (arbitrary units). |
sfp |
Sampling rate (Hz). |
a respiration object: time, volume (chest expansion, zero-mean),
flow (its time derivative), and rate.
rp <- respiration(duration = 30, rate = 15) rp$raterp <- respiration(duration = 30, rate = 15) rp$rate
Generates a heart-rate-variability tachogram by integral pulse frequency modulation: the instantaneous heart rate is modulated by respiration (faster on inspiration – respiratory sinus arrhythmia) and by a ~0.1 Hz Mayer wave, and a beat is emitted whenever the integral of the instantaneous rate crosses an integer. The resulting RR-interval series reproduces the high-frequency (respiratory) and low-frequency (Mayer) peaks of the HRV spectrum.
rsaTachogram( duration = 300, hr0 = 70, resp_rate = 15, rsa = 0.12, mayer = 0.06, mayer_freq = 0.1, ie_ratio = 0.5, noise = 0.01, seed = NULL )rsaTachogram( duration = 300, hr0 = 70, resp_rate = 15, rsa = 0.12, mayer = 0.06, mayer_freq = 0.1, ie_ratio = 0.5, noise = 0.01, seed = NULL )
duration |
Recording length (s). |
hr0 |
Mean heart rate (bpm). |
resp_rate |
Breathing rate (breaths/min); sets the HF peak. |
rsa |
Depth of the respiratory heart-rate modulation (fraction of |
mayer |
Depth of the Mayer-wave modulation (fraction of |
mayer_freq |
Mayer-wave frequency (Hz). |
ie_ratio |
Inspiration:expiration ratio of the driving respiration. |
noise |
SD of a white heart-rate jitter (fraction of |
seed |
Optional RNG seed. |
an rsaTachogram object: beat_time, rr (RR intervals, s), hr
(instantaneous, bpm), the driving resp waveform; and the HRV summaries
hf_peak (Hz of the HF spectral peak), lf_power, hf_power and lf_hf
(the LF/HF ratio).
Berger RD et al. (1989) IEEE Trans Biomed Eng 36:1049-1055 (IPFM). Task Force (1996) Circulation 93:1043-1065 (HRV bands).
hrv <- rsaTachogram(duration = 300, resp_rate = 15, seed = 1) hrv$hf_peak # ~ 0.25 Hz, the respiratory ratehrv <- rsaTachogram(duration = 300, resp_rate = 15, seed = 1) hrv$hf_peak # ~ 0.25 Hz, the respiratory rate
Runs each candidate intervention over the twin's PARAMETER POSTERIOR (sampled
from the twinReadiness() Fisher covariance), so every candidate gets a
predicted effect WITH a predictive interval and the probability that the effect
exceeds a minimal clinically important difference (MCID) – not a point estimate.
Candidates are ranked by a multi-criteria score combining efficacy, certainty
(probability of benefit), risk and burden, scaled by the twin's readiness
confidence; optional per-candidate evidence is reconciled in, and low-confidence
or evidence-disagreeing recommendations are flagged.
scoreInterventions( twin, candidates, readiness = NULL, outcome = "rom", mcid = NULL, n_mc = 200, weights = c(efficacy = 1, certainty = 1, risk = 1, burden = 0.5), duration = 6, dt = 0.005, theta0 = NULL, seed = NULL )scoreInterventions( twin, candidates, readiness = NULL, outcome = "rom", mcid = NULL, n_mc = 200, weights = c(efficacy = 1, certainty = 1, risk = 1, burden = 0.5), duration = 6, dt = 0.005, theta0 = NULL, seed = NULL )
twin |
A personalised |
candidates |
A list of |
readiness |
A |
outcome |
Outcome name ( |
mcid |
Minimal clinically important difference; defaults to 10% of the baseline outcome. |
n_mc |
Posterior samples for uncertainty propagation. |
weights |
Named weights |
duration, dt, theta0
|
Simulation settings; |
seed |
Optional RNG seed. |
an intervention_scores: a ranked ranking data frame (effect,
ci_lo/ci_hi, p_benefit, risk, burden, reconciled evidence, score,
flag), the mcid, and the readiness_confidence/verdict.
twinReadiness(), reconcileEvidence(), clinicalDecision()
truth <- limbTwin(strength = 6, damping = 0.6, stiffness = 3) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01) rdy <- twinReadiness(fit, imu, dt = 0.01) scoreInterventions(fit$twin, defaultInterventions(), rdy, n_mc = 50, seed = 1)$rankingtruth <- limbTwin(strength = 6, damping = 0.6, stiffness = 3) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01) rdy <- twinReadiness(fit, imu, dt = 0.01) scoreInterventions(fit$twin, defaultInterventions(), rdy, n_mc = 50, seed = 1)$ranking
Runs the movement twin and the cardiorespiratory generators and returns a small cross-modal summary: movement range of motion and peak velocity, and the mean heart rate.
simulateMultimodal( mt, duration_move = 6, duration_hr = 120, dt = 0.005, theta0 = 0.1, seed = NULL )simulateMultimodal( mt, duration_move = 6, duration_hr = 120, dt = 0.005, theta0 = 0.1, seed = NULL )
mt |
a |
duration_move |
Movement simulation length (s). |
duration_hr |
Heart-rate simulation length (s). |
dt, theta0
|
Movement simulation settings. |
seed |
Optional RNG seed. |
a named numeric vector: rom, peak_velocity, mean_hr, resp_rate.
multimodalTwin(), insilicoInterventionMultimodal()
simulateMultimodal(multimodalTwin(fitness = 0.6), seed = 1)simulateMultimodal(multimodalTwin(fitness = 0.6), seed = 1)
Integrates the twin's equation of motion with a fixed-step 4th-order Runge-Kutta scheme, returning the true joint state trajectory (angle, angular velocity and acceleration).
simulateTwin(twin, duration = 5, dt = 0.005, theta0 = 0, omega0 = 0)simulateTwin(twin, duration = 5, dt = 0.005, theta0 = 0, omega0 = 0)
twin |
A |
duration |
Simulation time (s). |
dt |
Time step (s). |
theta0, omega0
|
Initial angle and angular velocity. |
a twin_trajectory list: time, theta, omega, alpha (angular
acceleration), and twin.
simulateTwin(limbTwin(), duration = 2, dt = 0.005)$theta[1:5]simulateTwin(limbTwin(), duration = 2, dt = 0.005)$theta[1:5]
Apportions the variance of a model output among the parameters using the
Saltelli/Jansen estimator: first-order indices Si (each parameter's own
contribution) and total-effect indices STi (including interactions). A near-
zero total index means the parameter is non-influential (and typically
unidentifiable).
sobolIndices( ranges, output = NULL, base = limbTwin(), n = 256L, duration = 6, dt = 0.02, theta0 = 0.2, seed = NULL )sobolIndices( ranges, output = NULL, base = limbTwin(), n = 256L, duration = 6, dt = 0.02, theta0 = 0.2, seed = NULL )
ranges |
A named list of length-2 parameter ranges (uniform sampling). |
output |
A function |
base |
A |
n |
Base sample size (total model evaluations |
duration, dt, theta0
|
Simulation settings for the default output. |
seed |
Optional RNG seed. |
a sobol_indices data frame: param, first-order Si, total-effect
STi.
Saltelli A, et al. (2010) Comput Phys Commun 181:259-270; Jansen MJW (1999).
validateTwin(), profileLikelihood()
sobolIndices(list(strength = c(3, 9), damping = c(0.2, 1.0), stiffness = c(2, 6)), n = 64, seed = 1)sobolIndices(list(strength = c(3, 9), damping = c(0.2, 1.0), stiffness = c(2, 6)), n = 64, seed = 1)
The trust gate of the clinical loop. A personalised twin may only inform
decisions once it is shown to (1) FIT the patient's data and (2) have
IDENTIFIABLE parameters. This recomputes the personalisation residual (fit
quality) and a Fisher-information parameter covariance (identifiability), and
returns a verdict (ready / provisional / not_ready) with a [0, 1]
confidence that downstream scoring uses to scale and flag its recommendations.
twinReadiness( fit, imu, dt = 0.005, gyro_var = 0.03^2, accel_var = 0.15^2, cv_warn = 0.25, chi2_ready = 3, chi2_bad = 10 )twinReadiness( fit, imu, dt = 0.005, gyro_var = 0.03^2, accel_var = 0.15^2, cv_warn = 0.25, chi2_ready = 3, chi2_bad = 10 )
fit |
A |
imu |
The IMU data the twin was fitted to. |
dt |
Sample period (s). |
gyro_var, accel_var
|
Measurement-noise variances (match the fit). |
cv_warn |
Log-space parameter SD above which a parameter is deemed poorly identified. |
chi2_ready, chi2_bad
|
Reduced-chi-square thresholds: below |
a twin_readiness: verdict, reduced_chi2 (goodness of fit across
ALL channels), gyro_nrmse, the (misfit-inflated) parameter cov, per-
parameter cv, confidence, and the estimates/theta0/gyro_bias.
Sargent RG (2013) J Simul 7:12-24 (model validation); Bevington & Robinson (reduced chi-square goodness of fit); Cramer-Rao identifiability.
personalizeTwin(), scoreInterventions(), validateTwin()
truth <- limbTwin(strength = 7, damping = 0.5) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01) twinReadiness(fit, imu, dt = 0.01)$verdicttruth <- limbTwin(strength = 7, damping = 0.5) imu <- imuMeasure(simulateTwin(truth, 6, 0.01), imuSensor(gyro_noise = 0.03), seed = 1) fit <- personalizeTwin(imu, limbTwin(strength = 4), estimate = "strength", dt = 0.01) twinReadiness(fit, imu, dt = 0.01)$verdict
A sigma-point (Julier-Uhlmann) Kalman filter for a nonlinear state-space model
x_k = f(x_{k-1}, k) + w, y_k = h(x_k, k) + v. Used here to assimilate a
mechanistic twin with sensor data and to estimate augmented parameters, but it
is a stand-alone engine for any such model.
unscentedKalmanFilter( f, h, y, Q, R, x0, P0, alpha = 0.001, beta = 2, kappa = 0 )unscentedKalmanFilter( f, h, y, Q, R, x0, P0, alpha = 0.001, beta = 2, kappa = 0 )
f |
Process model |
h |
Measurement model |
y |
Measurement matrix ( |
Q |
Process-noise covariance (state x state). |
R |
Measurement-noise covariance (m x m). |
x0, P0
|
Initial state mean and covariance. |
alpha, beta, kappa
|
Sigma-point spread parameters (defaults 1e-3, 2, 0). |
a ukf_result: state (n_steps x L filtered means), cov (final
covariance), innovation (n_steps x m).
Julier SJ, Uhlmann JK (1997) SPIE 3068:182-193; Wan & van der Merwe (2000).
# track a noisy sine with a random-walk state set.seed(1); y <- matrix(sin(seq(0, 6, length.out = 60)) + rnorm(60, 0, .1), 60, 1) fit <- unscentedKalmanFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1)) fit$state[60]# track a noisy sine with a random-walk state set.seed(1); y <- matrix(sin(seq(0, 6, length.out = 60)) + rnorm(60, 0, .1), 60, 1) fit <- unscentedKalmanFilter(function(x, k) x, function(x, k) x, y, Q = matrix(0.05), R = matrix(0.01), x0 = 0, P0 = matrix(1)) fit$state[60]
Runs the applicable simulation-validity checks and returns a structured report: verification (energy conservation, integrator order), identifiability (can parameters be recovered in principle), Monte-Carlo recovery (bias / SD / coverage), sensitivity screening (which parameters matter), cross-condition prediction + face validity, and the domain of applicability. This is verification and internal validation on the model itself; validation against a subject's REAL data is done separately by personalising and predicting held-out recordings.
validateTwin( twin, estimate = c("strength", "damping"), sensor = imuSensor(gyro_noise = 0.03, accel_noise = 0.12), duration = 6, dt = 0.02, theta0 = 0.2, n_mc = 15L, n_traj = 8L, method = "optim", checks = c("verification", "identifiability", "recovery", "sensitivity", "prediction", "domain") )validateTwin( twin, estimate = c("strength", "damping"), sensor = imuSensor(gyro_noise = 0.03, accel_noise = 0.12), duration = 6, dt = 0.02, theta0 = 0.2, n_mc = 15L, n_traj = 8L, method = "optim", checks = c("verification", "identifiability", "recovery", "sensitivity", "prediction", "domain") )
twin |
A |
estimate |
Parameters whose identifiability/recovery are assessed. |
sensor |
The |
duration, dt, theta0
|
Simulation settings for the checks. |
n_mc, n_traj
|
Monte-Carlo recovery repetitions and Morris trajectories. |
method |
Assimilation method for the recovery check ( |
checks |
Which layers to run (subset of |
a twin_validity report.
Sargent RG (2013) J Simulation 7:12-24; ASME V&V 10/40; Oberkampf & Roy (2010) Verification and Validation in Scientific Computing.
personalizeTwin(), insilicoIntervention()
rep <- validateTwin(limbTwin(strength = 6, damping = 0.5, stiffness = 6), estimate = c("strength", "damping"), n_mc = 8, n_traj = 5) rep$verificationrep <- validateTwin(limbTwin(strength = 6, damping = 0.5, stiffness = 6), estimate = c("strength", "damping"), n_mc = 8, n_traj = 5) rep$verification
Simulates arterial pressure with the classic Windkessel model: a half-sine
ventricular ejection each beat charges the arterial compliance C, which
drains through the peripheral resistance R; the 3-element variant adds a
proximal characteristic impedance Zc. During diastole (no inflow) the
pressure decays exponentially with time constant R * C. With physiological
parameters it reproduces the systolic/diastolic/mean/pulse pressures and the
diastolic run-off.
windkessel( duration = 10, hr = 70, stroke_volume = 80, R = 1, C = 1.9, Zc = 0.05, systolic_frac = 0.3, elements = 3, sfp = 250, p0 = 80 )windkessel( duration = 10, hr = 70, stroke_volume = 80, R = 1, C = 1.9, Zc = 0.05, systolic_frac = 0.3, elements = 3, sfp = 250, p0 = 80 )
duration |
Recording length (s). |
hr |
Heart rate (bpm). |
stroke_volume |
Stroke volume ejected each beat (mL). |
R |
Peripheral (systemic vascular) resistance (mmHg*s/mL). |
C |
Arterial compliance (mL/mmHg). |
Zc |
Proximal characteristic impedance (mmHg*s/mL); used when |
systolic_frac |
Fraction of the cardiac cycle occupied by ejection. |
elements |
2 or 3 (2-element = pure RC; 3-element adds |
sfp |
Output sampling rate (Hz). |
p0 |
Initial compliance-node pressure (mmHg). |
a windkessel object: time, pressure (arterial pressure, mmHg),
flow (aortic inflow, mL/s), and the derived systolic, diastolic, map
(mean arterial pressure) and pp (pulse pressure), measured over the last
settled beats.
Westerhof N, Lankhaar J-W, Westerhof BE (2009) The arterial Windkessel. Med Biol Eng Comput 47:131-141.
baroreflex(), ecgsyn(), respiration()
bp <- windkessel(duration = 10, hr = 70) round(c(bp$systolic, bp$diastolic)) # ~ 120 / 80 mmHgbp <- windkessel(duration = 10, hr = 70) round(c(bp$systolic, bp$diastolic)) # ~ 120 / 80 mmHg