Healthcare

beta = 0.3; gam = 0.1;                 % infection and recovery rates
f = @(tt, y) [-beta*y(1)*y(2); beta*y(1)*y(2) - gam*y(2); gam*y(2)];
[tt, y] = ode45(f, [0 160], [0.99; 0.01; 0]);
plot(tt, y, 'LineWidth', 2)
legend('susceptible', 'infected', 'recovered')
xlabel('days'), ylabel('fraction of population')
title(sprintf('SIR model, R₀ = %.1f', beta/gam))

% A synthetic ECG with noise and baseline wander, and automatic R-peak detection
rng(3)
fs = 250; tt = 0:1/fs:8;
hr = 72; rr = 60/hr;                       % beats per minute, seconds per beat
ecg = zeros(size(tt));
for beat = 0:rr:8
  d = tt - beat;
  ecg = ecg + 0.12*exp(-((d - 0.10)/0.025).^2) ...   % P wave
            - 0.15*exp(-((d - 0.19)/0.008).^2) ...   % Q
            + 1.20*exp(-((d - 0.21)/0.010).^2) ...   % R
            - 0.25*exp(-((d - 0.23)/0.010).^2) ...   % S
            + 0.30*exp(-((d - 0.45)/0.045).^2);      % T wave
end
ecg = ecg + 0.04*randn(size(tt)) + 0.1*sin(2*pi*0.3*tt);

% R peaks: local maxima above a threshold
mid = ecg(2:end-1);
isPeak = [false, mid > ecg(1:end-2) & mid >= ecg(3:end) & mid > 0.6, false];
pk = find(isPeak);

plot(tt, ecg, 'LineWidth', 1), hold on
scatter(tt(pk), ecg(pk), 40, 'r', 'filled'), hold off
xlim([0 4])
xlabel('time (s)'), ylabel('mV')
title(sprintf('Detected heart rate: %.0f bpm', 60 / mean(diff(tt(pk)))))

% Survival in two simulated groups, followed for 36 months
function [x, y] = kmsteps(times, tmax)
  times = sort(times); n = numel(times);
  x = 0; y = 1;
  for i = 1:n
    if times(i) > tmax, break; end
    x = [x times(i) times(i)];
    y = [y y(end) 1 - i/n];
  end
  x = [x tmax]; y = [y y(end)];
end

rng(8)
control   = -log(rand(1, 150)) * 22;      % exponential survival times (months)
treatment = -log(rand(1, 150)) * 40;
[x1, y1] = kmsteps(control, 36);
[x2, y2] = kmsteps(treatment, 36);
plot(x1, y1, 'LineWidth', 2), hold on
plot(x2, y2, 'LineWidth', 2)
yline(0.5, 'k:', 'median')
hold off
ylim([0 1])
legend('control', 'treatment')
xlabel('months'), ylabel('survival probability')
title('Kaplan–Meier estimate')

% SIR model with a vaccinated fraction p of the population
beta = 0.32; gam = 0.1;                  % R0 = 3.2
for p = [0 0.3 0.5 0.65]
  f = @(tt, y) [-beta*y(1)*y(2); beta*y(1)*y(2) - gam*y(2)];
  [tt, y] = ode45(f, [0 200], [0.999 - p; 0.001]);
  plot(tt, 100*y(:, 2), 'LineWidth', 2, 'DisplayName', sprintf('%.0f%% vaccinated', 100*p))
  hold on
end
hold off
legend show
xlabel('days'), ylabel('infected (% of population)')
title('Herd immunity threshold: 1 − 1/R₀ ≈ 69%')