خوارزمية جيليسبي

في نظرية الاحتمالات ، تُولّد خوارزمية جيليسبي (أو خوارزمية دوب-جيليسبي أو خوارزمية المحاكاة العشوائية ، اختصارًا SSA ) مسارًا صحيحًا إحصائيًا (حلًا محتملًا) لنظام معادلات عشوائية تكون فيها معدلات التفاعل معروفة. ابتكرها جوزيف ل. دوب وآخرون (حوالي عام 1945)، وقدّمها دانيال جيليسبي عام 1976، وانتشر استخدامها عام 1977 في ورقة بحثية استخدمها فيها لمحاكاة أنظمة التفاعلات الكيميائية أو البيوكيميائية بكفاءة ودقة باستخدام قدرة حاسوبية محدودة (انظر المحاكاة العشوائية ). [ 1 ] مع تطور سرعة الحواسيب، استُخدمت الخوارزمية لمحاكاة أنظمة بالغة التعقيد. تُعدّ الخوارزمية مفيدة بشكل خاص لمحاكاة التفاعلات داخل الخلايا، حيث يكون عدد الكواشف منخفضًا، ويُمكن تتبّع كل تفاعل على حدة من الناحية الحسابية. رياضيًا، هي نوع من أنواع طريقة مونت كارلو الديناميكية ، وتُشبه طرق مونت كارلو الحركية . يُستخدم بكثرة في علم الأحياء الحاسوبي للأنظمة .

تاريخ

تتضمن العملية التي أدت إلى الخوارزمية عدة خطوات مهمة. ففي عام 1931، قدم أندريه كولموغوروف المعادلات التفاضلية التي تُعبر عن التطور الزمني للعمليات العشوائية التي تتم عبر القفزات، والمعروفة اليوم بمعادلات كولموغوروف (عملية ماركوف القفزية) (ويُعرف إصدار مبسط منها باسم المعادلة الرئيسية في العلوم الطبيعية). وفي عام 1940، اكتشف ويليام فيلر الشروط التي بموجبها تقبل معادلات كولموغوروف احتمالات (صحيحة) كحلول. وفي نظريته الأولى (1940)، أثبت أن الزمن اللازم للقفزة التالية يتوزع توزيعًا أُسّيًا، وأن احتمال وقوع الحدث التالي يتناسب مع معدل القفزة. وبذلك، أسس العلاقة بين معادلات كولموغوروف والعمليات العشوائية . وفي وقت لاحق، وسّع دوب (1942، 1945) حلول فيلر لتشمل حالات أخرى غير عمليات القفزات البحتة. طُبقت هذه الطريقة على الحواسيب بواسطة ديفيد جورج كيندال (1950) باستخدام حاسوب مانشستر مارك 1 ، ثم استخدمها موريس إس. بارتليت (1953) في دراساته حول تفشي الأوبئة. وقد توصل جيليسبي (1977) إلى الخوارزمية بطريقة مختلفة بالاعتماد على حجة فيزيائية.

فكرة

الرياضيات

يوجد في حجرة التفاعل عدد محدود من الجزيئات. وفي كل لحظة زمنية متناهية الصغر ، قد يحدث تفاعل واحد. ويتحدد معدل التفاعل بعدد الجزيئات في كل نوع كيميائي.

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

الافتراضات الرئيسية هي أن

  • كل رد فعل ماركوفي في الزمن
  • لا توجد علاقات بين ردود الفعل

بالنظر إلى الافتراضين، فإن وقت الانتظار العشوائي لبعض التفاعلات يتوزع بشكل أسي، حيث يكون المعدل الأسي هو مجموع معدلات التفاعلات الفردية.

الصلاحية في المحاكاة البيوكيميائية

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

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

الخوارزمية

تُقدّم مراجعة [ 1 ] ثلاثة صيغ مختلفة، ولكنها متكافئة؛ وهي: الطريقة المباشرة، وطريقة التفاعل الأول، وطريقة العائلة الأولى، حيث تُعتبر الطريقتان الأوليان حالتين خاصتين من الأخيرة. وترتكز صياغة الطريقتين المباشرة وطريقة التفاعل الأول على تنفيذ خطوات عكس مونت كارلو المعتادة على ما يُسمى "الفرضية الأساسية للحركية الكيميائية العشوائية "، والتي تُمثل رياضيًا الدالة

ص(τ،ج|x،ت)=أج(x)خبرة(-τجأج(x))،{\displaystyle p(\tau ,j\mid {\boldsymbol {x}},t)=a_{j}({\boldsymbol {x}})\exp \left(-\tau \sum _{j}a_{j}({\boldsymbol {x}})\right),}

حيث كل واحد منأ{\displaystyle a}المصطلحات هي دوال ميل لتفاعل أولي، وسيطها هوx{\displaystyle {\boldsymbol {x}}}، وهو متجه تعداد الأنواع.τ{\displaystyle \tau }المعامل هو الوقت اللازم للتفاعل التالي (أو وقت الإقامة)، وت{\displaystyle t}هو الوقت الحالي. وبتعبير آخر، يمكن قراءة هذا التعبير على أنه "الاحتمال، بالنظر إلى"X(ت)=x{\displaystyle {\boldsymbol {X}}(t)={\boldsymbol {x}}}، أي أن رد فعل النظام التالي سيحدث في فترة زمنية متناهية الصغر[ت+τ،ت+τ+دτ]{\displaystyle [t+\tau ,t+\tau +d\tau ]}وستكون ذات نسبة قياس العناصر المتوافقة معج{\displaystyle j}"التفاعل رقم 1". توفر هذه الصيغة نافذةً لطرق التفاعل المباشر والتفاعل الأول من خلال الإشارة إلىτ{\displaystyle \tau }هو متغير عشوائي ذو توزيع أسي ، وج{\displaystyle j}هو "متغير عشوائي صحيح مستقل إحصائياً باحتمالات نقطيةأج(x)/جأج(x){\displaystyle a_{j}({\boldsymbol {x}})/\sum _{j}a_{j}({\boldsymbol {x}})}".

وبالتالي، فإن طريقة توليد مونت كارلو هي ببساطة سحب رقمين شبه عشوائيين.ر1{\displaystyle r_{1}}ور2{\displaystyle r_{2}}على[0،1]{\displaystyle [0,1]}، وحساب

τ=1جأج(x)سجل(1ر1)،{\displaystyle \tau ={\frac {1}{\sum _{j}a_{j}({\boldsymbol {x}})}}\log \left({\frac {1}{r_{1}}}\right),}

و

ج={\displaystyle j={}}أصغر عدد صحيح يحققج=1جأج(x)>ر2جأج(x).{\displaystyle \sum _{j'=1}^{j}a_{j'}({\boldsymbol {x}})>r_{2}\sum _{j}a_{j}({\boldsymbol {x}}).}

باستخدام طريقة التوليد هذه لوقت الإقامة ورد الفعل التالي، تم صياغة خوارزمية الطريقة المباشرة بواسطة جيليسبي على النحو التالي:

1. تهيئة الوقتت=ت0{\displaystyle t=t_{0}}وحالة النظامx=x0{\displaystyle {\boldsymbol {x}}={\boldsymbol {x}}_{0}} 2. مع وجود النظام في حالةx{\displaystyle {\boldsymbol {x}}}في ذلك الوقتت{\displaystyle t}قم بتقييم جميعأج(x){\displaystyle a_{j}({\boldsymbol {x}})}ومجموعهاجأج(x){\displaystyle \sum _{j}a_{j}({\boldsymbol {x}})} 3. احسب القيمة المذكورة أعلاه لـτ{\displaystyle \tau }وج{\displaystyle j} 4. قم بإجراء التفاعل التالي عن طريق الاستبدالتت+τ{\displaystyle t\leftarrow t+\tau }وxx+νج{\displaystyle {\boldsymbol {x}}\leftarrow {\boldsymbol {x}}+\nu _{j}} 5. سجل(x،ت){\displaystyle ({\boldsymbol {x}},t)}حسب الرغبة. ارجع إلى الخطوة 2، وإلا أنهِ المحاكاة.

أينνج{\displaystyle \nu _{j}}يمثل إضافةجذ{\displaystyle j^{\text{th}}}مكون متجه تغيير الحالة المعطىν{\displaystyle \nu }تُعدّ هذه المجموعة من الخوارزميات مكلفة حسابيًا، ولذا توجد العديد من التعديلات والتحسينات عليها، بما في ذلك طريقة التفاعل التالي (جيبسون وبروك)، وتقنية قفز تاو ، بالإضافة إلى التقنيات الهجينة التي تُنمذج فيها المواد المتفاعلة الوفيرة بسلوك حتمي. تُؤثر التقنيات المُعدّلة عمومًا على دقة النظرية الكامنة وراء الخوارزمية عند ربطها بالمعادلة الرئيسية، لكنها تُقدّم تطبيقات معقولة لتحسين كبير في الأطر الزمنية. تُحدّد فئة اقتران شبكة التفاعلات التكلفة الحسابية للنسخ الدقيقة من الخوارزمية. في الشبكات ذات الاقتران الضعيف، يكون عدد التفاعلات المتأثرة بأي تفاعل آخر محدودًا بثابت صغير. أما في الشبكات ذات الاقتران القوي، فيمكن لتفاعل واحد أن يؤثر من حيث المبدأ على جميع التفاعلات الأخرى. وقد طُوّرت نسخة دقيقة من الخوارزمية ذات مقياس زمني ثابت للشبكات ذات الاقتران الضعيف، مما يُتيح محاكاة فعّالة للأنظمة ذات الأعداد الكبيرة جدًا من قنوات التفاعل (سليبوي طومسون بليمبتون 2008). طُوِّرت خوارزمية جيليسبي المعممة، التي تأخذ في الحسبان الخصائص غير الماركوڤية للأحداث البيوكيميائية العشوائية مع التأخير، بواسطة براتسون وآخرون (2005)، وبشكل مستقل بواسطة باريو وآخرون (2006)، وكذلك (كاي 2007). انظر المقالات المذكورة أدناه لمزيد من التفاصيل.

تتوفر صيغ الميل الجزئي، التي طُوّرت بشكل مستقل من قِبل كلٍّ من راماسوامي وآخرون (2009، 2010) وإندورخيا وبيال (2010)، لبناء مجموعة من النسخ الدقيقة للخوارزمية التي تتناسب تكلفتها الحسابية مع عدد الأنواع الكيميائية في الشبكة، بدلاً من عدد التفاعلات (الأكبر). يمكن لهذه الصيغ تقليل التكلفة الحسابية إلى مقياس زمني ثابت للشبكات ذات الترابط الضعيف، وإلى مقياس خطي على الأكثر مع عدد الأنواع للشبكات ذات الترابط القوي. كما اقتُرح أيضًا متغير ميل جزئي لخوارزمية جيليسبي المعممة للتفاعلات ذات التأخيرات (راماسوامي سبالزاريني 2011). يقتصر استخدام طرق الميل الجزئي على التفاعلات الكيميائية الأولية، أي التفاعلات التي تحتوي على متفاعلين مختلفين على الأكثر. يمكن تجزئة كل تفاعل كيميائي غير أولي بشكل مكافئ إلى مجموعة من التفاعلات الأولية، على حساب زيادة خطية (بحسب رتبة التفاعل) في حجم الشبكة.

أمثلة

ارتباط عكسي بين A و B لتكوين ثنائيات AB

قد يساعد مثال بسيط في شرح كيفية عمل خوارزمية جيليسبي. لنفترض نظامًا من جزيئات من نوعين، A و B. في هذا النظام، يرتبط A و B معًا بشكل عكسي لتكوين ثنائيات AB، بحيث يكون هناك تفاعلان ممكنان: إما أن يتفاعل A وB بشكل عكسي لتكوين ثنائي AB ، أو أن يتفكك ثنائي AB إلى A و B. ثابت معدل التفاعل لجزيء A واحد مع جزيء B واحد هوكد{\displaystyle k_{\mathrm {D} }}ومعدل تفاعل تفكك ثنائي AB هوكب{\displaystyle k_{\mathrm {B} }}.

إذا كان هناك جزيء واحد من كل نوع عند الزمن t ، فإن معدل تكوين ثنائي الوحدات هوكد{\displaystyle k_{\mathrm {D} }}أما إذا كان هناكنأ{\displaystyle n_{\mathrm {A} }}جزيئات من النوع أ ونب{\displaystyle n_{\mathrm {B} }}بالنسبة لجزيئات النوع B ، يكون معدل تكوين ثنائي الوحدات هوكدنأنب{\displaystyle k_{\mathrm {D}}n_{\mathrm {A} }n_{\mathrm {B} }}إذا كان هناكنأب{\displaystyle n_{\mathrm {AB} }}إذا كانت ثنائيات، فإن معدل تفكك الثنائيات هوكبنأب{\displaystyle k_{\mathrm {B} }n_{\mathrm {AB} }}.

معدل التفاعل الكلي،Rتيياتي{\displaystyle R_{\mathrm {TOT} }}، عند الزمن t يُعطى بالعلاقة التالية

Rتيياتي=كدنأنب+كبنأب{\displaystyle R_{\mathrm {TOT}}=k_{\mathrm {D}}n_{\mathrm {A} }n_{\mathrm {B} }+k_{\mathrm {B} }n_{\mathrm {AB} }}

لقد وصفنا الآن نموذجًا بسيطًا يتضمن تفاعلين. هذا التعريف مستقل عن خوارزمية جيلسبي. سنشرح الآن كيفية تطبيق خوارزمية جيلسبي على هذا النظام.

في هذه الخوارزمية، نتقدم زمنيًا على مرحلتين: حساب الوقت اللازم للتفاعل التالي، وتحديد أي من التفاعلات الممكنة هو التفاعل التالي. يُفترض أن التفاعلات عشوائية تمامًا، لذا إذا كان معدل التفاعل عند الزمن t هوRتيياتي{\displaystyle R_{\mathrm {TOT} }}ثم يكون الوقت، δt ، حتى حدوث التفاعل التالي هو رقم عشوائي يتم سحبه من دالة التوزيع الأسي بمتوسط1/Rتيياتي{\displaystyle 1/R_{\mathrm {TOT} }}وبالتالي، فإننا نتقدم بالوقت من t إلى t + δ t .

رسم بياني يوضح عدد جزيئات A (المنحنى الأسود) وثنائيات AB كدالة للزمن. بما أننا بدأنا بـ 10 جزيئات من A و B عند الزمن t = 0، فإن عدد جزيئات B يساوي دائمًا عدد جزيئات ولذلك لم يتم عرضه.

إن احتمال أن يكون هذا التفاعل عبارة عن ارتباط جزيء A بجزيء B هو ببساطة نسبة المعدل الكلي الناتج عن هذا النوع من التفاعل، أي

احتمالية حدوث رد الفعل هيP(أ+بAB)=كدنأنب/Rإجمالي{\displaystyle P({\ce {{A}+ B -> AB}})=k_{D}n_{A}n_{B}/R_{{\ce {TOT}}}}

إن احتمال أن يكون التفاعل التالي هو تفكك ثنائي AB هو 1 ناقص ذلك. لذا، باستخدام هذين الاحتمالين، إما أن نُشكّل ثنائيًا عن طريق الاختزالنأ{\displaystyle n_{\mathrm {A} }}ونب{\displaystyle n_{\mathrm {B} }}بزيادة واحد، وزيادةنأب{\displaystyle n_{\mathrm {AB} }}بمقدار واحد، أو نقوم بتفكيك ثنائي الوحدات وزيادةنأ{\displaystyle n_{\mathrm {A} }}ونب{\displaystyle n_{\mathrm {B} }}بمقدار واحد وينقصنأب{\displaystyle n_{\mathrm {AB} }}واحداً تلو الآخر.

الآن، قمنا بتمديد الزمن إلى t + δt ، وأجرينا تفاعلًا واحدًا. تقوم خوارزمية جيليسبي بتكرار هاتين الخطوتين عدة مرات حسب الحاجة لمحاكاة النظام للمدة التي نريدها (أي لعدد التفاعلات المطلوب). نتيجة محاكاة جيليسبي التي تبدأ بـنأ=نب=10{\displaystyle n_{\mathrm {A} }=n_{\mathrm {B} }=10}ونأب=0{\displaystyle n_{\mathrm {AB} }=0}عند t = 0، وحيثكد=2{\displaystyle k_{\mathrm {D} }=2}وكب=1{\displaystyle k_{\mathrm {B} }=1}يظهر ذلك على اليمين. بالنسبة لقيم هذه المعلمات، يوجد في المتوسط ​​8نأب{\displaystyle n_{\mathrm {AB} }}تُعدّ ثنائيات الجزيئات وجزيئتان من A و B ، ولكن نظرًا لقلة عدد الجزيئات، تكون التقلبات حول هذه القيم كبيرة. غالبًا ما تُستخدم خوارزمية جيليسبي لدراسة الأنظمة التي تكون فيها هذه التقلبات مهمة.

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

مراجع

  1. 1 2 جيليسبي، دانيال ت. (2007-05-01). "المحاكاة العشوائية للحركية الكيميائية" . المراجعة السنوية للكيمياء الفيزيائية . 58 (1): 35-55 . Bibcode : 2007ARPC...58...35G . doi : 10.1146/annurev.physchem.58.032806.104637 . ISSN 0066-426X . PMID 17037977 .  

للمزيد من القراءة