كيفية نمذجة الوباء باستخدام R عبر نموذج SEIR

مقدمة: لماذا تُعد نمذجة الأوبئة مهمة؟

أصبحت epidemiology أو علم الأوبئة من أكثر المجالات أهمية في العصر الحديث، لأنها تدرس كيفية تأثير الأمراض على التجمعات السكانية، وخاصة الأمراض المعدية مثل COVID-19. ومن أهم الأدوات المستخدمة لفهم انتشار العدوى ما يُعرف باسم نمذجة الأوبئة، وهي عملية بناء نماذج كمية تساعد على وصف مسار المرض والتنبؤ بتطوره بمرور الوقت.

المنهج الكلاسيكي في هذا المجال يعتمد على ما يسمى النماذج التقسيمية أو compartmental models. تقوم هذه الفكرة على تقسيم السكان إلى مجموعات وفق الحالة الصحية، ثم تحديد معدلات الانتقال بين هذه المجموعات، وبعد ذلك تُصاغ معادلات تفاضلية تصف تطور الوباء.

أنواع النماذج التقسيمية في علم الأوبئة

نموذج SI: أبسط تمثيل للعدوى

يُعد نموذج SI أبسط أشكال النمذجة الوبائية، إذ يقسم السكان إلى مجموعتين فقط:

• Susceptible: أفراد قابلون للإصابة.
• Infectious: أفراد مصابون وقادرون على نقل العدوى.

نموذج SIR: إضافة المتعافين

يضيف نموذج SIR قسماً ثالثاً هو Recovered، أي الأفراد الذين تعافوا من العدوى. ويُستخدم هذا النموذج كثيراً كنقطة انطلاق أساسية في الدراسات الوبائية لأنه يصف تطور الوباء بشكل مبسط لكنه مفيد.

نموذج SEIR: تمثيل أكثر واقعية

يُعد نموذج SEIR تطويراً مباشراً لـ SIR، إذ يضيف فئة جديدة هي Exposed لتمثيل الأشخاص الذين تعرضوا للفيروس لكنهم لم يصبحوا معديين بعد. وهذا يجعل النموذج أكثر دقة في وصف أمراض تمتلك فترة حضانة.

ما هو نموذج SEIR؟

يتكون النموذج الأساسي SEIR من أربع حجرات رئيسية:

• Susceptible: أفراد لم يتعرضوا للفيروس بعد.
• Exposed: أفراد تعرضوا للفيروس لكنهم لم يصبحوا ناقلين للعدوى حتى الآن.
• Infectious: أفراد أصبحوا قادرين على نقل المرض.
• Recovered: أفراد تعافوا واكتسبوا مناعة.

ويُرمز إلى حجم المجتمع الكلي بالرمز N، وهو مجموع الأفراد في جميع هذه الأقسام.

المعاملات الأساسية في نموذج SEIR

معامل الانتقال β أو beta

يمثل β معامل انتقال العدوى، ويمكن النظر إليه على أنه متوسط عدد المخالطات المعدية التي يُحدثها الفرد المصاب خلال وحدة زمنية. كلما ارتفعت قيمة β زادت فرص انتشار الفيروس.

معدل التحول إلى العدوى σ أو sigma

يمثل σ المعدل الذي يتحول به الأفراد من حالة Exposed إلى حالة Infectious. وهو يساوي مقلوب متوسط الزمن اللازم ليصبح الشخص معدياً. فإذا كان المتوسط 4 أيام، فإن σ = 1/4 أي 0.25.

معدل التعافي γ أو gamma

يعبر γ عن المعدل الذي يتعافى به الأفراد المصابون. وإذا كان متوسط وقت التعافي 10 أيام، فإن γ = 1/10 أي 0.1.

معدل الوفيات μ أو mu

يُستخدم μ كمعامل اختياري لتمثيل معدل الوفاة بين المصابين. وكلما ارتفعت قيمته زادت فتك العدوى ضمن النموذج.

المعادلات التفاضلية في نموذج SEIR

بعد تحديد الأقسام والمعاملات، يمكن بناء نظام من المعادلات التفاضلية يصف تغير كل فئة بمرور الزمن.

المعادلة الأولى: الفئة المعرضة للإصابة S

لا يمكن أن تزداد قيمة S في هذا النموذج، لأنها لا تستقبل تدفقاً من أقسام أخرى. بل تتناقص فقط عندما ينتقل الأفراد إلى حالة التعرض للعدوى. وبافتراض امتزاج مثالي للسكان، تكون صيغة التغير كما يلي:

المعادلة الثانية: الفئة المتعرضة E

تستقبل فئة E الأفراد الخارجين من S، ثم يفقد هذا القسم أفراده عندما يصبحون معديين بمعدل σ.

المعادلة الثالثة: الفئة المعدية I

ينتقل الأفراد إلى القسم I من القسم E. ثم يغادرون هذا القسم بطريقتين: التعافي بمعدل γ أو الوفاة بمعدل μ.

المعادلة الرابعة: الفئة المتعافية R

تزداد قيمة R مع تعافي الأفراد المصابين، ولا تغادرها الحالات في الصيغة الأساسية للنموذج.

المعادلة الخامسة: الوفيات M بشكل اختياري

يمكن أيضاً إضافة قسم خامس يمثل الوفيات. وإذا ضبطنا μ على القيمة 0 فيمكن تجاهل هذا الجزء بالكامل.

وبهذا نحصل على النظام الكامل للمعادلات:

العدد التكاثري الأساسي R₀

من أهم المؤشرات في أي نموذج وبائي ما يسمى العدد التكاثري الأساسي أو R₀، وهو مقياس يقدّر عدد الأشخاص الذين قد ينقل إليهم الفرد المصاب العدوى في المتوسط.

• إذا كانت قيمة R₀ أكبر من 1، فمن المرجح أن يتحول التفشي إلى وباء واسع الانتشار.
• إذا كانت قيمة R₀ أقل من 1، فغالباً يمكن احتواء العدوى.

كيفية حل معادلات SEIR باستخدام لغة R

للاستفادة من النموذج عملياً، لا بد من حل المعادلات بالنسبة للزمن. ويمكن تنفيذ ذلك باستخدام لغة R ومكتبة deSolve المتخصصة في حل المعادلات التفاضلية.

تعريف دالة النموذج

require(deSolve)

SEIR <- function(time, current_state, params){
  with(as.list(c(current_state, params)),{
    N <- S + E + I + R
    dS <- -(beta * S * I) / N
    dE <- (beta * S * I) / N - sigma * E
    dI <- sigma * E - gamma * I - mu * I
    dR <- gamma * I
    dM <- mu * I
    return(list(c(dS, dE, dI, dR, dM)))
  })
}

يستورد هذا المقطع الحزمة deSolve، ثم يعرّف دالة باسم SEIR(). تستقبل الدالة ثلاثة مدخلات:

1. الزمن الحالي time.
2. الحالة الحالية للنظام current_state.
3. قائمة المعاملات params.

وفي داخل الدالة يتم تعريف المعادلات التفاضلية ثم إرجاعها بنفس ترتيب المتغيرات.

تهيئة المعاملات والحالة الابتدائية

params <- c(beta = 0.5, sigma = 0.25, gamma = 0.2, mu = 0.001)
initial_state <- c(S = 999999, E = 1, I = 0, R = 0, M = 0)
times <- 0:365

في هذا المثال نبدأ بعدد سكان يساوي مليون شخص تقريباً، مع وجود حالة تعرض واحدة فقط في البداية، ثم نتابع تطور الوباء خلال 365 يوماً.

حل النموذج باستخدام ode()

model <- ode(initial_state, times, SEIR, params)

تستخدم الدالة ode() من حزمة deSolve لحل المعادلات مع الزمن اعتماداً على:

• الحالة الابتدائية لكل قسم.
• متجه الزمن.
• الدالة SEIR().
• المعاملات اللازمة.

عرض ملخص النتائج

summary(model)

        S            E               I               R            M
Min.    108263.6   3.616607e-07    0.000000e+00    0.00        0.0000
1st Qu. 108263.7   5.957435e-03    1.414971e-02    63894.43    319.4721
Median  108395.7   8.470071e+00    1.273726e+01    886814.36   4434.0718
Mean    362798.6   9.745754e+03    1.212158e+04    612272.74   3061.3637
3rd Qu. 852375.5   1.734331e+03    2.533956e+03    887299.83   4436.4991
Max.    999999.0   1.092967e+05    1.265161e+05    887299.86   4436.4993
N       366.0      3.660000e+02    3.660000e+02    366.00      366.0000
sd      381257.2   2.475783e+04    2.969234e+04    387333.47   1936.6673

تُظهر هذه النتائج عدداً من الملاحظات المهمة:

• نحو 108,264 شخصاً من أصل مليون لم يُصابوا بالعدوى.
• بلغت الذروة حوالي 126,516 حالة معدية في الوقت نفسه.
• وصل عدد المتعافين في نهاية الفترة إلى نحو 887,300.
• بلغ إجمالي الوفيات حوالي 4,436 حالة.

رسم تطور الوباء

matplot(model, type = "l", lty = 1, main = "SEIR model", xlab = "Time")
legend <- colnames(model)[2:6]
legend("right", legend = legend, col = 2:6, lty = 1)

يمكن كذلك استخدام مكتبات متقدمة مثل ggplot2 للحصول على رسوم بيانية أكثر احترافية.

تحليل ذروة الإصابات

infections <- as.data.frame(model)$I
peak <- max(infections)
match(peak, infections)

يكشف هذا التحليل أن ذروة الإصابات حدثت في اليوم 112 من بداية النموذج.

نمذجة التدخلات الحكومية مثل الإغلاق

يفترض نموذج SEIR الأساسي أن سلوك المجتمع لا يتغير مع الوقت. لكن الواقع مختلف، إذ تتدخل الحكومات أو يتغير سلوك الأفراد بشكل مباشر لتقليل المخالطات. ويمكن محاكاة هذا الأثر عبر تعديل قيمة β.

تعريف نموذج SEIR مع الإغلاق

SEIR_lockdown <- function(time, current_state, params){
  with(as.list(c(current_state, params)),{
    beta = ifelse((time <= start_lockdown || time >= end_lockdown), 0.5, 0.1)
    N <- S + E + I + R
    dS <- -(beta * S * I) / N
    dE <- (beta * S * I) / N - sigma * E
    dI <- sigma * E - gamma * I - mu * I
    dR <- gamma * I
    dM <- mu * I
    return(list(c(dS, dE, dI, dR, dM)))
  })
}

الاختلاف الوحيد هنا هو استخدام الدالة ifelse() لخفض قيمة β إلى 0.1 خلال فترة الإغلاق، ما يعكس انخفاض فرص انتقال العدوى.

إعداد فترة الإغلاق وتشغيل النموذج

params <- c(sigma = 0.25, gamma = 0.2, mu = 0.001, start_lockdown = 90, end_lockdown = 150)
initial_state <- c(S = 999999, E = 1, I = 0, R = 0, M = 0)
times <- 0:365
model <- ode(initial_state, times, SEIR_lockdown, params)

في هذا السيناريو يبدأ الإغلاق في اليوم 90 وينتهي في اليوم 150.

ملخص نتائج النموذج مع التدخل

summary(model)

        S            E               I            R            M
Min.    156885.7   7.699207e-01    0.00000      0.00        0.0000
1st Qu. 160478.2   6.929205e+01    97.71405     63668.75    318.3438
Median  789214.4   1.246389e+03    1735.66330   194379.16   971.8958
Mean    589558.9   9.216918e+03    11460.62036  387824.44   1939.1222
3rd Qu. 867639.6   1.030043e+04    13780.17591  829898.56   4149.4928
Max.    999999.0   6.083432e+04    72443.97892  838916.89   4194.5845
N       366.0      3.660000e+02    366.00000    366.00      366.0000
sd      350719.3   1.570278e+04    18893.31145  346542.57   1732.7128

وتشير النتائج إلى ما يلي:

• نحو 156,886 شخصاً لم يُصابوا بالعدوى.
• بلغت الذروة حوالي 72,444 حالة معدية متزامنة.
• وصل عدد المتعافين إلى 838,917 شخصاً.
• بلغ إجمالي الوفيات نحو 4,195 حالة.

يتضح هنا أن التدخل خفّض الذروة مقارنة بالنموذج الأساسي، لكنه قد يؤدي أيضاً إلى ظهور موجة ثانية بعد انتهاء الإغلاق.

رسم النموذج مع التدخل

matplot(model, type = "l", lty = 1, main = "SEIR model (with intervention)", xlab = "Time")
legend <- colnames(model)[2:6]
legend("right", legend = legend, col = 2:6, lty = 1)

تحديد يوم الذروة الجديدة

infections <- as.data.frame(model)$I
peak <- max(infections)
match(peak, infections)

في هذا السيناريو حدثت ذروة الإصابات في اليوم 223.

كيف يمكن تطوير النموذج أكثر؟

يمكن توسيع النموذج الأساسي بعدة طرق للحصول على تمثيل أقرب إلى الواقع، مثل:

• إضافة أثر اللقاحات على انتقال العدوى.
• تغيير قيمة β موسمياً.
• تقسيم السكان إلى فئات عمرية.
• تمثيل الاختلافات الجغرافية بين المدن والمناطق.
• إدخال معدلات ولادة ووفيات طبيعية في النماذج طويلة الأمد.

قيود نموذج SEIR والنماذج التقسيمية

على الرغم من فائدة نموذج SEIR، فإنه يعتمد على افتراضات تبسيطية قد لا تصمد دائماً أمام الواقع المعقد. من أبرز هذه القيود:

• يفترض امتزاجاً متجانساً بين السكان، بينما المخالطة الفعلية غالباً محلية وغير متساوية.
• يفترض مجتمعاً معزولاً، في حين أن التنقل بين المناطق قد يعيد إدخال الفيروس مراراً.
• لا يتعامل في صورته الأساسية مع البنية العمرية للسكان.
• يعتمد على متوسطات للمعاملات، بينما الواقع يشهد تفاوتاً كبيراً بين الأفراد.
• لا يصف بدقة الظواهر الدقيقة على النطاق الصغير، مثل أحداث العدوى الفائقة super-spreading.

لذلك، فإن هذا النموذج مناسب لفهم الاتجاه العام للوباء على مستوى كلي، لكنه ليس أداة مثالية لتوقع التفاصيل الدقيقة جداً.

لماذا تبقى هذه النماذج مفيدة رغم القيود؟

تكمن قوة هذه النماذج في بساطتها وقدرتها على تقديم إطار رياضي واضح لفهم حركة المرض. كما أن توفر أدوات مثل R وPython ومكتباتهما الإحصائية يجعل من السهل على الباحثين والمحللين البدء في استكشاف هذه التقنيات وتطويرها.

الخلاصة التقنية

يُعد نموذج SEIR من أفضل النماذج التمهيدية لفهم ديناميكيات انتشار الأوبئة، لأنه يوازن بين البساطة والقدرة التفسيرية. واستخدام لغة R مع مكتبة deSolve يمنح المحلل وسيلة عملية لحل المعادلات ورسم السيناريوهات المختلفة، بما في ذلك أثر التدخلات مثل الإغلاق أو خفض المخالطة. تقنياً، لا ينبغي التعامل مع هذه النماذج على أنها تنبؤات قطعية، بل كأدوات تحليل تساعد على استيعاب السلوك العام للعدوى واتخاذ قرارات مبنية على تصور كمي أفضل.