I. LỊCH SỬ Vào giữa thế kỷ 14, ở châu Âu xuất hiện đại dịch “Cái Chết Đen” quét qua lục địa này. Những nghiên cứu gần đây cho thấy “Cái Chết Đen” hình thành do loài bọ chét truyền từ chuột sang người, gây nên bệnh dịch hạch. Đại dịch này xuất phát ở Ý vào khoảng tháng 12, năm 1347 từ những chuyến tàu cập bến từ phương Đông. Vài năm sau đó, đại dịch lan rộng khoảng 200 đến 400 dặm mỗi năm, khiến cho một phần ba dân số châu Âu thiệt mạng và 80% những người nhiễm bệnh sẽ chết trong 2 đến 3 ngày tới. Hình dưới đây mô tả mức độ lan truyền theo dạng sóng của bệnh dịch này. Ảnh minh hoạ sự lan truyền đại dịch “Cái Chết Đen” ở châu Âu vào năm 1347 – 1350 Trong bài này, ta nghiên cứu một mô hình đơn giản của sự lan truyền và cách xấp xỉ trong thực tế với giả định rằng lượng tổng thể sinh vật luôn cố định. Giả sử có một dịch bệnh lan truyền trong một khu vực, bệnh này có thể được cứu chữa, nhưng người mắc bệnh có thể chết, ta không loại bỏ số lượng người chết ra khỏi tổng dân số trong vùng. Ta phân tổng số người trong vùng đó thành 3 lớp: - Lớp dễ bệnh $S$: Những người trong lớp này chưa hề mắc bệnh và có nguy cơ nhiễm bệnh. - Lớp nhiễm bệnh $I$: Những người trong lớp này đã mắc bệnh và có khả năng truyền bệnh sang người khác. - Lớp hết bệnh $R$: Những người trong lớp này đã được trị khỏi bệnh hoặc đã chết vì bệnh. Đây là mô hình $SIR$ mà ta sẽ nghiên cứu sau đây. II. MÔ HÌNH SIR 2.1. Điều kiện nghiên cứu - Bệnh dịch xảy ra trong khoảng thời gian đủ ngắn để lượng dân số luôn cố định, tính cả người đã chết vì bệnh dịch này vào tổng số dân. - Chu kỳ ủ bệnh không đáng kể. - Nếu người nhiễm bệnh đã hết bệnh thì người này không còn khả năng nhiễm bệnh. Do đó, người này sẽ ở lại lớp $R$. - Mẫu dân số đủ lớn để có kết quả xấp xỉ đúng. - Ta xác định mức độ lan truyền dịch bệnh bằng định luật tác dụng khối lượng như sau: $$S+I \overset{r}{\mathop{\to }} \,2I$$ $$I \overset{a}{\mathop{\to }} \,R$$ Ý nghĩa: $S+I\overset{r}{\mathop{\to }}\,2I$: Ở vế trái, người trong lớp $S$ bị người trong lớp $I$ lây bệnh với tốc độ $r>0$, khiến người đó chuyển sang lớp $I$, thu được vế phải là $I+I=2I$. $I\overset{a}{\mathop{\to }}\,R$: Người trong lớp $I$ sau một thời gian sẽ hết bệnh (hoặc chết vì bệnh này) và chuyển sang lớp $R$ với tốc độ $a>0$. 2.2. Mô hình Ta xem mỗi lớp là một hàm số theo thời gian $t$ gồm $S\left( t \right),I\left( t \right)$ và $R\left( t \right)$ có tính chất: (i) Lớp nhiễm bệnh có tốc độ tỉ lệ thuận với số lượng người nhiễm bệnh và người dễ bệnh, tức $rSI$, với $r>0$ là tham số hằng, đó cũng là tốc độ mất đi số người trong lớp dễ bệnh. (ii) Tốc độ hết bệnh của người nhiễm bệnh tỉ lệ thuận với số lượng người nhiễm bệnh, tức $aI$, với $a>0$ là hằng số, $1/a$ là độ đo thời gian một người ở trong trạng thái nhiễm bệnh. (iii) Chu kỳ ủ bệnh ngắn, tức người dễ bệnh khi tiếp xúc với mầm bệnh sẽ nhiễm bệnh ngay. Từ các tính chất trên, ta có được mô hình cổ điển Kermack – McKendrick, do W. O. Kermack và A. G. McKendrick công bố vào năm 1927. $$\frac{dS}{dt}=-rSI \, \, \, (1)$$ $$\frac{dI}{dt}=rSI-aI \, \, \, (2)$$ $$\frac{dR}{dt}=aI \, \, \, (3)$$ với $r>0$ là tốc độ lây nhiễm và $a>0$ là tốc độ hết bệnh của người nhiễm bệnh. Ta muốn tìm nghiệm không âm với mỗi $S,I$ và $R$, tức $S\left( t \right),I\left( t \right),R\left( t \right)\geq 0$ với mọi $t\geq 0$. Do tổng kích thước dân số luôn cố định nên ta cộng các phương trình từ (1) đến (3) thu được: $$\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}=0$$ $$\Rightarrow S\left( t \right)+I\left( t \right)+R\left( t \right)=N \, \, \, (4)$$ với $N$ là tổng kích thước dân số. Giả sử ta có điều kiện đầu $S\left( t=0 \right)={{S}_{0}},~I\left( t=0 \right)={{I}_{0}},~R\left( t=0 \right)=0$, tức tại thời điểm bắt đầu khảo sát, ta có ${{S}_{0}}$ người chưa mắc bệnh, ${{I}_{0}}$ người đã nhiễm bệnh và chưa có ai khỏi bệnh (hoặc chết vì bệnh đó). Từ phương trình (4), ta được: $$\frac{d}{dt}\left( S+I+R \right)=0$$ $$S\left( 0 \right)+I\left( 0 \right)+R\left( 0 \right)={{S}_{0}}+{{I}_{0}}$$ Trong bất kỳ dịch bệnh nào, người ta thường hỏi rằng nếu biết $r,a,{{S}_{0}}$ và ${{I}_{0}}$ người nhiễm bệnh thì 1. Dịch bệnh có lan rộng nữa không? Mức độ lan rộng theo thời gian như thế nào? Khi nào dịch bệnh kết thúc? Từ phương trình (1) $$\frac{dS}{dt}=-rSI$$ Do vế phải không thể dương nên $S$ giảm, hay nói cách khác, số người dễ bệnh giảm theo thời gian. Do đó $S\leq {{S}_{0}}$. Mặt khác, từ phương trình (2) $$\frac{dI}{dt}=rSI-aI=I\left( rS-a \right)<I\left( r{{S}_{0}}-a \right)$$ Vậy dịch bệnh sẽ không còn lan rộng, ít nhất là giữ nguyên nếu $I\left( r{{S}_{0}}-a \right)<0$ và dịch bệnh sẽ lan rộng nếu $I\left( r{{S}_{0}}-a \right)>0$. Nếu ${{S}_{0}}<a/r$ thì $$\frac{dI}{dt}=I\left( rS-a \right)\leq 0, \forall t\geq 0$$ Khi đó ${{I}_{0}}>I\left( t \right)\to 0$ khi $t\to \infty $, tức số lượng người nhiễm bệnh ngày càng giảm, dẫn đến dịch bệnh không còn khả năng lan rộng. Mặt khác, nếu ${{S}_{0}}>a/r$ thì $I\left( t \right)$ tăng, dẫn đến dịch bệnh lan rộng. Khái niệm “lan rộng” có nghĩa rằng $I\left( t \right)>{{I}_{0}}$ với một vài giá trị $t>0$, từ đó ta có khái niệm ngưỡng hiện tượng, nếu ${{S}_{0}}>{{S}_{c}}=a/r$ thì dịch bệnh lan rộng, còn nếu ${{S}_{0}}<{{S}_{c}}$ thì không. Tham số tới hạn $\rho =a/r$ còn được gọi là tốc độ hết bệnh tương đối. 2. Nếu dịch bệnh lan rộng thì có nhiều nhất bao nhiêu người sẽ nhiễm bệnh ứng với thời gian cho trước bất kỳ? Từ biểu thức định luật tác dụng khối lượng $S+I\overset{r}{\mathop{\to}}\,2I$, ta nhận thấy sự thay đổi số người nhiễm bệnh có phụ thuộc đến sự thay đổi số người dễ bệnh. Do đó, ta lấy phương trình (2) chia cho phương trình (1), thu được tốc độ thay đổi tức thời của người nhiễm bệnh theo người dễ bệnh $$\frac{dI}{dS}=-\frac{\left( rS-a \right)I}{rSI}=-1+\frac{\rho }{S};$$ $$\rho =\frac{a}{r},~\left( I\ne 0 \right) \, \, \, \, (5)$$ Lấy tích phân phương trình (5), ta được mặt phẳng quỹ đạo pha $\left( I,~S \right)$, hình dưới đây là đồ thị của mặt phẳng này: Nghiệm số của mô hình SIR của phương trình (1) - (3), đường nét đứt là $S+I={{S}_{0}}+{{I}_{0}}$, đường trơn là quỹ đạo pha, $r=0.01,a=0.25$ $$\frac{dI}{dS}=-1+\frac{\rho }{S}$$ $$\Leftrightarrow dI=\left( -1+\frac{\rho }{S} \right)dS$$ $$\Rightarrow \int dI=\int \left( -1+\frac{\rho }{S} \right)dS$$ $$\Rightarrow I+ \text{ const }=-S+\rho \ln S+\text{ const }$$ $$\Rightarrow I+S-\rho \ln S=\text{ const }={{I}_{0}}+{{S}_{0}}-\rho \ln {{S}_{0}} \, \, \, \, (6)$$ Để tìm giá trị cực đại, trước hết ta tìm nghiệm của phương trình $dI/dS=0$ $$\frac{dI}{dS}=0\Leftrightarrow -1+\frac{\rho }{S}=0\Rightarrow S=\rho$$ Trong trường hợp ${{S}_{0}}<\rho $, khi đó dịch bệnh không lan rộng, số lượng người nhiễm bệnh có xu hướng giảm xuống. Từ đó, ta được ${{I}_{\text{max}}}={{I}_{0}}$, tức số người nhiễm bệnh nhiều nhất xác định tại $t=0$ với ${{I}_{0}}$ người. Trong trường hợp ${{S}_{0}}>\rho $, thay $S=\rho $ vào phương trình (6), ta được: $$I+\rho -\rho \ln \rho ={{I}_{0}}+{{S}_{0}}-\rho \ln {{S}_{0}}$$ $$\Leftrightarrow I={{I}_{0}}+{{S}_{0}}-\rho \ln {{S}_{0}}+\rho \ln \rho -\rho$$ $$\Leftrightarrow I=N-\rho +\rho \ln \left( \frac{\rho }{{{S}_{0}}} \right) \, \, \, \, (7)$$ Biểu thức (7) của $I$ chính là giá trị lớn nhất, do đó ta có kết luận về số người nhiễm bệnh nhiều nhất như sau: $${{I}_{\text{max}}}=\left\{ \begin{matrix} {{I}_{0}} \, \, \, ;{{S}_{0}}<\rho \\ {{I}_{0}}+{{S}_{0}}-\rho \ln {{S}_{0}}-\rho \ln \rho -\rho \, \, \,; {{S}_{0}}>\rho \\ \end{matrix} \right.$$ 3. Tổng cộng có bao nhiêu người nhiễm bệnh? Vì trục $I=0$ là đường kỳ dị nên mọi quỹ đạo $I$ tiến đến $0$ khi $t\to\infty$. Mặt khác, người trong lớp $I$ sau một thời gian sẽ hết bệnh (hoặc chết vì bệnh này) và chuyển sang lớp $R$. Do đó, tổng số người mắc bệnh bằng $R\left( \infty \right)$ xác định bằng cách chuyển vế phương trình (4) với $t=\infty $ như sau $$R\left( \infty \right)=N-S\left( \infty \right)-I\left( \infty \right)=N-S\left( \infty \right) \, \, \, \, (8)$$ Ta tính $S\left( \infty \right)$ bằng cách thay $t=\infty $ vào phương trình (6) như sau $$I\left( \infty \right)+S\left( \infty \right)-\rho \ln S\left( \infty \right)={{I}_{0}}+{{S}_{0}}-\rho \ln {{S}_{0}}$$ $$\Leftrightarrow S\left( \infty \right)-\rho \ln S\left( \infty \right)=N-\rho \ln {{S}_{0}}$$ Ta giải phương trình (8) để tìm $S\left( \infty \right)$, với $0\leq S\left( \infty \right)\leq {{S}_{0}}$. Phương trình (8) luôn có nghiệm (chứng minh bằng cách đặt $x=S\left( \infty \right)$, với $0\leq x\leq {{S}_{0}}$), thay nghiệm này vào phương trình (8), ta thu được tổng số người nhiễm bệnh. 2.3. Ứng dụng của mô hình SIR Vào năm 1978, tạp chí y khoa Anh The Lancet có bài viết cung cấp một số dữ liệu thống kê về dịch cúm tại một trường nội trú nam sinh. Trường này có tổng cộng 763 nam sinh, từ ngày 22/1/1978 đến ngày 4/2/1978 có tổng cộng 512 em nhiễm bệnh và người ta dự đoán ban đầu có 1 em bị bệnh. Từ dữ kiện trên, ta được $N=763,{{S}_{0}}=762,I=1.\rho =202,r=2.18\times {{10}^{-3}}/\text{ngày}$, xấp xỉ phương trình (1)-(3) được nghiệm $S\left( t \right)$ và $I\left( t \right)$, nghiệm $R\left( t \right)$ tỷ lệ với diện tích dưới đường cong $I\left( t \right)$. Đồ thị nghiệm xấp xỉ của $S\left( t \right)$ và $I\left( t \right)$ với dữ liệu dịch cúm (dấm chấm tròn đen) từ tạp chí The Lancet Mô hình SIR còn được áp dụng trong không gian không thuần nhất, nơi mà vùng diện tích chịu ảnh hưởng bởi dịch bệnh luôn thay đổi. Để mô hình hoá và ước lượng nghiệm xấp xỉ, người ta sử dụng đến sóng lưu động.
