Mô hình toán học về bệnh truyền nhiễm đóng vai trò quan trọng trong y tế công cộng. Trong đại dịch COVID-19, các mô hình phương trình vi phân được sử dụng rộng rãi để dự báo xu hướng và đánh giá biện pháp kiểm soát. Bài viết này hướng dẫn chi tiết cách xây dựng mô hình SEID (Susceptible-Exposed-Infectious-Recovered) bằng Matlab, từ phương trình đến trực quan hóa kết quả.

1. Lý thuyết mô hình SEID và chuẩn bị môi trường

Mô hình SEID mở rộng từ SIR, thêm trạng thái người nhiễm tiềm ẩn (E). Nhiều bệnh như cúm, sởi có giai đoạn ủ bệnh khi người nhiễm chưa có triệu chứng nhưng đã có khả năng lây nhiễm. Mô hình chia dân số thành bốn nhóm:

• S (Susceptible): Người khỏe, có nguy cơ nhiễm bệnh
• E (Exposed): Người đã nhiễm nhưng chưa phát bệnh
• I (Infectious): Người bệnh có triệu chứng, có thể lây truyền
• R (Recovered): Người đã hồi phục và có miễn dịch

Hệ phương trình vi phân cơ bản:

dS/dt = -β×S×I/N
dE/dt = β×S×I/N - σ×E
dI/dt = σ×E - γ×I
dR/dt = γ×I

Trong đó N = S+E+I+R là tổng dân số, β là tỷ lệ lây nhiễm, σ là tỷ lệ chuyển từ tiềm ẩn sang nhiễm, γ là tỷ lệ hồi phục.

1.1 Kiểm tra môi trường Matlab

Yêu cầu: Matlab R2018a trở lên, Toolbox tối ưu hóa (Optimization Toolbox).

% Kiểm tra toolbox
if ~license('test','Optimization_Toolbox')
    error('Cần cài đặt Optimization Toolbox');
end

2. Giải hệ phương trình vi phân bằng ode45

ode45 là solver phổ biến nhất trong Matlab, sử dụng phương pháp Runge-Kutta bậc 4-5 với bước thích ứng. Cần chuyển hệ phương trình SEID thành hàm Matlab.

2.1 Định nghĩa hàm phương trình vi phân

Tạo file seid_system.m:

function dy = seid_system(t, y, p)
    beta = p.beta;
    sigma = p.sigma;
    gamma = p.gamma;
    
    S = y(1);
    E = y(2);
    I = y(3);
    R = y(4);
    N = sum(y(1:4));
    
    dS = -beta * S * I / N;
    dE = beta * S * I / N - sigma * E;
    dI = sigma * E - gamma * I;
    dR = gamma * I;
    
    dy = [dS; dE; dI; dR];
end

2.2 Thiết lập điều kiện ban đầu và giải

Giả sử thành phố 1 triệu dân, ban đầu có 10 người tiềm ẩn và 5 người nhiễm:

% Điều kiện ban đầu [S0, E0, I0, R0]
y0 = [1e6-15, 10, 5, 0];  

% Thời gian mô phỏng (ngày)
t_span = [0 200];  

% Tham số
p.beta = 0.3;      % Tỷ lệ tiếp xúc/ngày
p.sigma = 1/5;     % Thời gian ủ bệnh 5 ngày
p.gamma = 1/10;    % Thời gian nhiễm 10 ngày

% Giải
[t, y] = ode45(@(t,y) seid_system(t,y,p), t_span, y0);

2.3 Trực quan hóa kết quả

figure('Position', [100 100 800 500])
plot(t, y(:,1), 'b-', 'LineWidth', 2); hold on;
plot(t, y(:,2), 'm--', 'LineWidth', 2);
plot(t, y(:,3), 'r-.', 'LineWidth', 2);
plot(t, y(:,4), 'g:', 'LineWidth', 2);
xlabel('Thời gian (ngày)'); ylabel('Số người');
legend('S (Dễ nhiễm)','E (Tiềm ẩn)','I (Nhiễm)','R (Hồi phục)');
title('Diễn biến mô hình SEID');
grid on;
set(gca, 'FontSize', 12);

3. Phân tích độ nhạy và hiệu chỉnh tham số

Trong thực tế, tham số mô hình cần được hiệu chỉnh dựa trên dữ liệu quan sát. Sử dụng phương pháp bình phương tối thiểu.

3.1 Chuẩn bị dữ liệu quan sát

% Dữ liệu quan sát [thời gian, số ca nhiễm]
obs = [
    0    5
    30   25
    60   180
    90   650
    120  320
    150  120
    180  50
];

3.2 Hàm mục tiêu

function err = calc_error(p, t_span, y0, obs)
    [t, y] = ode45(@(t,y) seid_system(t,y,p), t_span, y0);
    sim_I = interp1(t, y(:,3), obs(:,1));
    err = sum((sim_I - obs(:,2)).^2);
end

3.3 Tối ưu hóa tham số

% Giá trị khởi tạo
guess = struct('beta', 0.2, 'sigma', 0.25, 'gamma', 0.1);

% Tùy chọn tối ưu
opts = optimset('Display', 'iter', 'MaxIter', 100);

% Chạy tối ưu
best_p = fminsearch(@(p) calc_error(p, t_span, y0, obs), guess, opts);

disp('Tham số tối ưu:');
disp(best_p);

4. Mở rộng mô hình

Mô hình cơ bản có thể mở rộng theo đặc điểm bệnh cụ thể. Dưới đây là hai biến thể phổ biến.

4.1 Thêm tiêm chủng

Giả sử mỗi ngày có tỷ lệ v người dễ nhiễm được tiêm chủng:

function dy = seid_vaccine(t, y, p)
    beta = p.beta;
    sigma = p.sigma;
    gamma = p.gamma;
    v = p.v;  % Tỷ lệ tiêm chủng/ngày
    
    S = y(1);
    E = y(2);
    I = y(3);
    R = y(4);
    N = sum(y(1:4));
    
    dS = -beta*S*I/N - v*S;
    dE = beta*S*I/N - sigma*E;
    dI = sigma*E - gamma*I;
    dR = gamma*I + v*S;
    
    dy = [dS; dE; dI; dR];
end

4.2 Tỷ lệ tiếp xúc thay đổi theo thời gian

function dy = seid_time_varying(t, y, p)
    beta = p.beta0 * exp(-p.k * t);  % Giảm dần theo thời gian
    sigma = p.sigma;
    gamma = p.gamma;
    
    S = y(1);
    E = y(2);
    I = y(3);
    R = y(4);
    N = sum(y(1:4));
    
    dS = -beta*S*I/N;
    dE = beta*S*I/N - sigma*E;
    dI = sigma*E - gamma*I;
    dR = gamma*I;
    
    dy = [dS; dE; dI; dR];
end

5. Code hoàn chỉnh và mẹo gỡ lỗi

Script tích hợp tất cả chức năng trên:

% Mô hình SEID hoàn chỉnh
clear; close all; clc;

% Tham số
p.beta = 0.35;
p.sigma = 1/5;
p.gamma = 1/10;
p.v = 0.001;  % Tùy chọn tiêm chủng

% Điều kiện ban đầu
N = 1e6;
y0 = [N-15, 10, 5, 0];

% Thời gian
t_span = [0 200];

% Giải
[t, y] = ode45(@(t,y) seid_vaccine(t,y,p), t_span, y0);

% Vẽ
figure('Position', [100 100 900 600])
subplot(2,1,1)
plot(t, y(:,1), 'b-', t, y(:,2), 'm--', t, y(:,3), 'r-.', t, y(:,4), 'g:')
legend('S','E','I','R'); xlabel('Ngày'); ylabel('Số người');
title('Diễn biến mô hình SEID');

% Số ca mới hàng ngày
new_cases = p.sigma * y(:,2);
subplot(2,1,2)
plot(t, new_cases, 'k-', 'LineWidth', 2)
xlabel('Ngày'); ylabel('Số ca mới/ngày');
title('Dự báo ca nhiễm mới hàng ngày');

Mẹo gỡ lỗi thường gặp:

• Lỗi giải phương trình: Kiểm tra hàm trả về vector cột, tham số trong khoảng hợp lý (β ∈ [0,1]).
• Kết quả không hợp lý: Giảm bước thời gian: options = odeset('MaxStep', 0.1);
• Tối ưu không hội tụ: Chuẩn hóa tham số về [0,1], thử nhiều giá trị khởi tạo khác nhau.